library(pracma)
library(vegan)## Loading required package: permute
## Loading required package: lattice
## This is vegan 2.5-3
##
## Attaching package: 'vegan'
## The following object is masked from 'package:pracma':
##
## procrustes
library(dplyr)##
## Attaching package: 'dplyr'
## The following objects are masked from 'package:stats':
##
## filter, lag
## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, union
library(data.table)##
## Attaching package: 'data.table'
## The following objects are masked from 'package:dplyr':
##
## between, first, last
library(ggplot2)
library(RColorBrewer)
library("affy") ## Loading required package: BiocGenerics
## Loading required package: parallel
##
## Attaching package: 'BiocGenerics'
## The following objects are masked from 'package:parallel':
##
## clusterApply, clusterApplyLB, clusterCall, clusterEvalQ,
## clusterExport, clusterMap, parApply, parCapply, parLapply,
## parLapplyLB, parRapply, parSapply, parSapplyLB
## The following objects are masked from 'package:dplyr':
##
## combine, intersect, setdiff, union
## The following objects are masked from 'package:stats':
##
## IQR, mad, sd, var, xtabs
## The following objects are masked from 'package:base':
##
## anyDuplicated, append, as.data.frame, cbind, colMeans,
## colnames, colSums, do.call, duplicated, eval, evalq, Filter,
## Find, get, grep, grepl, intersect, is.unsorted, lapply,
## lengths, Map, mapply, match, mget, order, paste, pmax,
## pmax.int, pmin, pmin.int, Position, rank, rbind, Reduce,
## rowMeans, rownames, rowSums, sapply, setdiff, sort, table,
## tapply, union, unique, unsplit, which, which.max, which.min
## Loading required package: Biobase
## Welcome to Bioconductor
##
## Vignettes contain introductory material; view with
## 'browseVignettes()'. To cite Bioconductor, see
## 'citation("Biobase")', and for packages 'citation("pkgname")'.
library("limma")##
## Attaching package: 'limma'
## The following object is masked from 'package:BiocGenerics':
##
## plotMA
library("multtest")
library("DESeq2")## Loading required package: S4Vectors
## Loading required package: stats4
##
## Attaching package: 'S4Vectors'
## The following objects are masked from 'package:data.table':
##
## first, second
## The following objects are masked from 'package:dplyr':
##
## first, rename
## The following object is masked from 'package:base':
##
## expand.grid
## Loading required package: IRanges
##
## Attaching package: 'IRanges'
## The following object is masked from 'package:data.table':
##
## shift
## The following objects are masked from 'package:dplyr':
##
## collapse, desc, slice
## Loading required package: GenomicRanges
## Loading required package: GenomeInfoDb
## Loading required package: SummarizedExperiment
## Loading required package: DelayedArray
## Loading required package: matrixStats
##
## Attaching package: 'matrixStats'
## The following objects are masked from 'package:Biobase':
##
## anyMissing, rowMedians
## The following object is masked from 'package:dplyr':
##
## count
##
## Attaching package: 'DelayedArray'
## The following objects are masked from 'package:matrixStats':
##
## colMaxs, colMins, colRanges, rowMaxs, rowMins, rowRanges
## The following object is masked from 'package:base':
##
## apply
library(scales)
library(tidyr)##
## Attaching package: 'tidyr'
## The following object is masked from 'package:S4Vectors':
##
## expand
library("hgu95av2cdf")## Warning: replacing previous import 'AnnotationDbi::tail' by 'utils::tail'
## when loading 'hgu95av2cdf'
## Warning: replacing previous import 'AnnotationDbi::head' by 'utils::head'
## when loading 'hgu95av2cdf'
##
library("pheatmap")
library("dendextend")##
## ---------------------
## Welcome to dendextend version 1.12.0
## Type citation('dendextend') for how to cite the package.
##
## Type browseVignettes(package = 'dendextend') for the package vignette.
## The github page is: https://github.com/talgalili/dendextend/
##
## Suggestions and bug-reports can be submitted at: https://github.com/talgalili/dendextend/issues
## Or contact: <tal.galili@gmail.com>
##
## To suppress this message use: suppressPackageStartupMessages(library(dendextend))
## ---------------------
##
## Attaching package: 'dendextend'
## The following object is masked from 'package:data.table':
##
## set
## The following object is masked from 'package:permute':
##
## shuffle
## The following object is masked from 'package:stats':
##
## cutree
library(funrar)
library(edgeR)
setwd("~/Desktop/spinal cord injury project/Ion torrent paper")
Colvec_expedition=read.csv("Color_vector.csv",header=T,row.names=1,sep=",")L2=read.csv("otu_table.L2.csv",header=T,row.names=1,sep=",")
L2=t(L2)
L2=t(L2)
#calculate the relative abundances of bacterial strains
L2_abundances=make_relative(L2)
colvec_all <- c("#999999", "#999999","#999999","#999999","#999999","#E69F00","#E69F00","#E69F00","#E69F00","#E69F00", "#56B4E9","#56B4E9","#56B4E9","#56B4E9","#56B4E9")
# Applying a log transformation
d.clr_L2= log2(L2_abundances+1)
# Calculating Bray-Curtis dissimilarity
Bray_L2=vegdist(d.clr_L2,"bray")
# ANOSIM statistics
L2_anosim=anosim(Bray_L2, permutations = 999, grouping=Colvec_expedition$Zones)
summary(L2_anosim)##
## Call:
## anosim(x = Bray_L2, grouping = Colvec_expedition$Zones, permutations = 999)
## Dissimilarity: bray
##
## ANOSIM statistic R: 0.3973
## Significance: 0.002
##
## Permutation: free
## Number of permutations: 999
##
## Upper quantiles of permutations (null model):
## 90% 95% 97.5% 99%
## 0.134 0.186 0.252 0.303
##
## Dissimilarity ranks between and within classes:
## 0% 25% 50% 75% 100% N
## Between 4 37.50 59.0 80.50 105 75
## Lam 1 3.50 7.0 12.50 20 10
## T10 9 15.75 23.5 44.25 84 10
## T4 7 59.75 81.0 90.25 98 10
#ANOSIM statistic R: 0.3973 , Significance: 0.001
# Applying a non-constrained ordination (PCoA) on the dissimilarity matrix
L2_bray=capscale(Bray_L2~-1)
# Extracting the percentage explained by the first two dimensions and automatically adding them to the axes titles
L2_bray_eig = eigenvals(L2_bray)
percentage_variance_explained <- L2_bray_eig / sum(L2_bray_eig)
sum_percentage_variance_explained <- cumsum(L2_bray_eig / sum(L2_bray_eig))
xlabel= as.numeric(format(round((percentage_variance_explained[1]*100), 2), nsmall = 2))
xlabel= sprintf("%.2f %%", xlabel)
xlabel= paste ("PCo1 (", xlabel, "of variation explained", ")")
ylabel= as.numeric(format(round((percentage_variance_explained[2]*100), 2), nsmall = 2))
ylabel= sprintf("%.2f %%", ylabel)
ylabel= paste ("PCo2 (", ylabel,"of variation explained", ")")
# Plotting the figure
# Adding the axes, grid, and other aestethics
plot(L2_bray, type="n", xlab="", xlim=c(-1,1),ylab="", tck = -0.01, mgp = c(3, 0.2, 0),
xaxp = c(-4, 4, 8), panel.first=grid(col = "white",lty=0))
title(ylab=ylabel, line=2, cex.lab=1)
title(xlab=xlabel, line=2, cex.lab=1)
abline(h=0, v=0, col = "white", lty = 1, lwd = 1.5)
abline(h=-10:10, v=-10:10, col = "lightgray", lty = "dotted", lwd = 0.5)
par(lty=2)
points(L2_bray,cex= 2.5,pch=21, col="black", bg= colvec_all, lwd = 1)
# Adding a legend
legend(-2,1, pt.cex=2.5 , pt.lwd = 1, c("Lam", "T10", "T4"), bty = "n", pch = 21, col="black", pt.bg = c("#999999","#E69F00","#56B4E9"), cex = 1.5)
ordihull(L2_bray,Colvec_expedition$Zones, display = "sites", col = NULL,label=T)Bray_L2_SCI=read.csv("Bray_L2_SCI.csv",header=T,sep=",")
pairwise.wilcox.test(Bray_L2_SCI$Bray_curtis, Bray_L2_SCI$Zones,p.adjust.method = "fdr")##
## Pairwise comparisons using Wilcoxon rank sum test
##
## data: Bray_L2_SCI$Bray_curtis and Bray_L2_SCI$Zones
##
## T10
## T4 0.0052
##
## P value adjustment method: fdr
#p-value = 0.0052
boxplot_Bray_L2_SCI=ggplot(Bray_L2_SCI, aes(x=Zones, y=Bray_curtis)) + geom_boxplot(fill=c( "#E69F00", "#56B4E9"))+geom_jitter(width = 0.06)+labs(title=NULL,x=NULL,y =c( "Bray_curtis dissimilarity"))+theme_classic()+geom_point(shape=16,fill="black",size=2)
boxplot_Bray_L2_SCI1=boxplot_Bray_L2_SCI+theme(plot.title = element_text(size = 28,hjust=0.5))+theme(axis.text.x = element_text(color="black",size = 20,face = "bold"))+theme(axis.text.y = element_text(color="black",size = 18),axis.title.y=element_text(size=20))
#Including boxplots
p=boxplot_Bray_L2_SCI1+scale_y_continuous(expand = c(0,0), limits=c(0.3,0.9))+geom_segment(y=0.83,yend=0.83,x=1,xend=2)+geom_text(x=1.5,y=0.83,label="**",size=10)+theme(plot.margin=unit(c(1,1,1.5,1.2),"cm"))+theme(axis.title.y = element_text(margin = margin(t = 0, r = 20, b = 0, l = 0)))
p#Read bacterial abundances in phylum level
phylum_singleM=read.csv("15_ion_annatable_phylum_singleM.csv",header=T,row.names=1,sep=",")
phylum_singleM1=t(phylum_singleM)
phylum_singleM2=as.data.frame(phylum_singleM1)
phylum_singleM2$Sample=factor(rownames(phylum_singleM2))
phylum_singleM2_m=melt(phylum_singleM2)## Using Sample as id variables
phylum_singleM_p=ggplot(phylum_singleM2_m)+geom_bar(aes(x=Sample,y=value,fill=variable),stat="identity")+labs(title = NULL,x=NULL,y = "Relative abundance")+ scale_y_continuous(limits = c(0,1.00), expand = c(0, 0))+theme_bw()
# legend labels
phylum_singleM_p=phylum_singleM_p+geom_segment(y=0,yend=1.00,x=5.5,xend=5.5)+geom_segment(y=0,yend=1.00,x=10.5,xend=10.5)
phylum_singleM_p+theme(plot.title = element_text(size = 20,hjust=0.5))+scale_x_discrete(labels=c("1","2","3","4","5","6","7","8","9","10","11","12","13","14","15"))+theme(axis.text.x = element_text(color="black",size = 12))+guides(fill=guide_legend(title=NULL))+theme(legend.title = element_text(size=18))+ theme(legend.text = element_text(size=14),legend.position="right")+guides(fill=guide_legend("Phylum"))+scale_fill_brewer(palette = "Set2")#On the x-axis, the numbers 1–15 represent individual mice within each group. numbers 1-5 represents Lam, numbers 6-10 represents T10 and numbers 11-15 represents T4##Wilcoxon rank sum test for bacterial phyla
clin.data <- data.frame(ID = 1:ncol(phylum_singleM), group = rep(c("Lam", "T10","T4"), each=5))
clin.data## ID group
## 1 1 Lam
## 2 2 Lam
## 3 3 Lam
## 4 4 Lam
## 5 5 Lam
## 6 6 T10
## 7 7 T10
## 8 8 T10
## 9 9 T10
## 10 10 T10
## 11 11 T4
## 12 12 T4
## 13 13 T4
## 14 14 T4
## 15 15 T4
idx_Lam<-sample(which(clin.data $group=="Lam"), 5, replace=FALSE)
idx_T10 <- sample(which(clin.data $group=="T10"), 5, replace=FALSE)
idx_T4 <- sample(which(clin.data $group=="T4"), 5, replace=FALSE)
#Wilcoxon rank sum test between Lam and T4 at phylum level
phylum_data_test <- phylum_singleM[,c(idx_Lam, idx_T4)]
clin.data_phylum_sub <- data.frame(ID = 1:ncol(phylum_data_test), group = rep(c("Lam","T4"), each=5))
phylum_data_test1=t(phylum_data_test)
phylum_data_test2=cbind(phylum_data_test1,clin.data_phylum_sub)
wilcox.test(phylum_data_test2$` p__Bacteroidetes`~phylum_data_test2$group)##
## Wilcoxon rank sum test
##
## data: phylum_data_test2$` p__Bacteroidetes` by phylum_data_test2$group
## W = 19, p-value = 0.2222
## alternative hypothesis: true location shift is not equal to 0
#p-value = 0.22
wilcox.test(phylum_data_test2$` p__Firmicutes`~phylum_data_test2$group)##
## Wilcoxon rank sum test
##
## data: phylum_data_test2$` p__Firmicutes` by phylum_data_test2$group
## W = 14, p-value = 0.8413
## alternative hypothesis: true location shift is not equal to 0
#p-value = 0.8413
wilcox.test(phylum_data_test2$` p__Firmicutes_A`~phylum_data_test2$group)##
## Wilcoxon rank sum test
##
## data: phylum_data_test2$` p__Firmicutes_A` by phylum_data_test2$group
## W = 5, p-value = 0.1508
## alternative hypothesis: true location shift is not equal to 0
#p-value = 0.1508
wilcox.test(phylum_data_test2$` p__Actinobacteria`~phylum_data_test2$group)## Warning in wilcox.test.default(x = c(0.004, 0.008, 0.007, 0.005, 0.012), :
## cannot compute exact p-value with ties
##
## Wilcoxon rank sum test with continuity correction
##
## data: phylum_data_test2$` p__Actinobacteria` by phylum_data_test2$group
## W = 2.5, p-value = 0.04653
## alternative hypothesis: true location shift is not equal to 0
#p-value = 0.04653 *
wilcox.test(phylum_data_test2$` p__Proteobacteria`~phylum_data_test2$group)## Warning in wilcox.test.default(x = c(0, 0.001, 0.001, 0.001, 0.002), y =
## c(0.001, : cannot compute exact p-value with ties
##
## Wilcoxon rank sum test with continuity correction
##
## data: phylum_data_test2$` p__Proteobacteria` by phylum_data_test2$group
## W = 12.5, p-value = 1
## alternative hypothesis: true location shift is not equal to 0
#p-value = 1
##Wilcoxon rank sum test obetween Lam and T10 at phylum level
phylum_data_test4 <- phylum_singleM[,c(idx_Lam, idx_T10)]
clin.data_phylum_sub1 <- data.frame(ID = 1:ncol(phylum_data_test), group = rep(c("Lam","T10"), each=5))
phylum_data_test5=t(phylum_data_test4)
phylum_data_test6=cbind(phylum_data_test5,clin.data_phylum_sub1)
wilcox.test(phylum_data_test6$` p__Bacteroidetes`~phylum_data_test6$group)##
## Wilcoxon rank sum test
##
## data: phylum_data_test6$` p__Bacteroidetes` by phylum_data_test6$group
## W = 13, p-value = 1
## alternative hypothesis: true location shift is not equal to 0
#p-value = 1
wilcox.test(phylum_data_test6$` p__Firmicutes`~phylum_data_test6$group)##
## Wilcoxon rank sum test
##
## data: phylum_data_test6$` p__Firmicutes` by phylum_data_test6$group
## W = 23, p-value = 0.03175
## alternative hypothesis: true location shift is not equal to 0
#p-value = 0.03175 *
wilcox.test(phylum_data_test6$` p__Firmicutes_A`~phylum_data_test6$group)##
## Wilcoxon rank sum test
##
## data: phylum_data_test6$` p__Firmicutes_A` by phylum_data_test6$group
## W = 11, p-value = 0.8413
## alternative hypothesis: true location shift is not equal to 0
#p-value = 0.8413
wilcox.test(phylum_data_test6$` p__Actinobacteria`~phylum_data_test6$group)## Warning in wilcox.test.default(x = c(0.004, 0.008, 0.007, 0.005, 0.012), :
## cannot compute exact p-value with ties
##
## Wilcoxon rank sum test with continuity correction
##
## data: phylum_data_test6$` p__Actinobacteria` by phylum_data_test6$group
## W = 1.5, p-value = 0.0278
## alternative hypothesis: true location shift is not equal to 0
#p-value = 0.0278 *
wilcox.test(phylum_data_test6$` p__Proteobacteria`~phylum_data_test6$group)## Warning in wilcox.test.default(x = c(0, 0.001, 0.001, 0.001, 0.002), y =
## c(0, : cannot compute exact p-value with ties
##
## Wilcoxon rank sum test with continuity correction
##
## data: phylum_data_test6$` p__Proteobacteria` by phylum_data_test6$group
## W = 20.5, p-value = 0.08326
## alternative hypothesis: true location shift is not equal to 0
# p-value = 0.08326#Wilcoxon rank sum test corrected by fdr method
singleM_phylum=cbind(phylum_singleM1,Colvec_expedition)
colnames(singleM_phylum)=c("Bacteroidetes","Firmicutes","Firmicutes_A","Actinobacteria","Proteobacteria","Firmicutes_B","others","Zones")
pairwise.wilcox.test(singleM_phylum$Actinobacteria, singleM_phylum$Zones,p.adjust.method = "fdr")## Warning in wilcox.test.default(xi, xj, paired = paired, ...): cannot
## compute exact p-value with ties
## Warning in wilcox.test.default(xi, xj, paired = paired, ...): cannot
## compute exact p-value with ties
## Warning in wilcox.test.default(xi, xj, paired = paired, ...): cannot
## compute exact p-value with ties
##
## Pairwise comparisons using Wilcoxon rank sum test
##
## data: singleM_phylum$Actinobacteria and singleM_phylum$Zones
##
## Lam T10
## T10 0.07 -
## T4 0.07 0.53
##
## P value adjustment method: fdr
pairwise.wilcox.test(singleM_phylum$Firmicutes, singleM_phylum$Zones,p.adjust.method = "fdr")##
## Pairwise comparisons using Wilcoxon rank sum test
##
## data: singleM_phylum$Firmicutes and singleM_phylum$Zones
##
## Lam T10
## T10 0.095 -
## T4 0.841 0.464
##
## P value adjustment method: fdr
#Boxplots for Firmicutes
Firmicutes=ggplot(singleM_phylum, aes(x=Zones, y=Firmicutes)) + geom_boxplot(fill=c("#999999", "#E69F00", "#56B4E9"))+geom_dotplot(binaxis='y', stackdir='center', dotsize=0.6)+labs(title = NULL,x=NULL,y = "Firmicutes \n relative abundances")+theme_classic()
Firmicutes_boxplot=Firmicutes+theme(plot.title = element_text(size = 20,hjust=0.5))+theme(axis.text.x = element_text(color="black",size = 20))+theme(axis.text.y = element_text(color="black",size = 16),axis.title.y=element_text(size=18))
Firmicutes_boxplot+scale_y_continuous(expand = c(0,0), limits = c(0,0.25))+geom_segment(y=0.23,yend=0.23,x=1,xend=2)+geom_text(x=1.5,y=0.23,label="*",size=12)+
theme(plot.margin=unit(c(1,1,1.5,1.2),"cm"))+theme(axis.title.y = element_text(margin = margin(t = 0, r = 20, b = 0, l = 0)))## `stat_bindot()` using `bins = 30`. Pick better value with `binwidth`.
Actinobacteria=ggplot(singleM_phylum, aes(x=Zones, y=Actinobacteria)) + geom_boxplot(fill=c("#999999", "#E69F00", "#56B4E9"))+geom_dotplot(binaxis='y', stackdir='center', dotsize=0.6)+labs(title = NULL,x=NULL,y = "Actinobacteria \n relative abundances")+theme_classic()
Actinobacteria_boxplot=Actinobacteria+theme(plot.title = element_text(size = 20,hjust=0.5))+theme(axis.text.x = element_text(color="black",size = 20))+theme(axis.text.y = element_text(color="black",size = 16),axis.title.y=element_text(size=18))
Actinobacteria_boxplot+scale_y_continuous(expand = c(0,0), limits = c(0,0.035))+geom_segment(y=0.032,yend=0.032,x=1,xend=2)+geom_text(x=1.5,y=0.032,label="*",size=12)+geom_segment(y=0.026,yend=0.026,x=1,xend=3)+geom_text(x=2,y=0.026,label="*",size=12)+theme(plot.margin=unit(c(1,1,1.5,1.2),"cm"))+theme(axis.title.y = element_text(margin = margin(t = 0, r = 20, b = 0, l = 0)))## `stat_bindot()` using `bins = 30`. Pick better value with `binwidth`.
#*p = 0.05 by Wilcoxon rank sum test . Each dot in the boxplot represents an individual mouse sample, and each group (Lam, T4, T10) contains five samples. #Read bacterial abundances in genus
genus_singleM=read.csv("15_ion_annatable_genus_singleM.csv",header=T,row.names=1,sep=",")
##Wilcoxon rank sum test
#Lam and T4
genus_singleM1<- genus_singleM[,c(idx_Lam, idx_T4)]
clin.sub_genus_singleM1 <- data.frame(ID = 1:ncol(genus_singleM1), group = rep(c("Lam","T4"), each=5))
fx <- function(x) {
res <-wilcox.test(x ~clin.sub_genus_singleM1$group)
return(res$p.value)
}
genus_singleM.wilcox1 <-apply(as.matrix(genus_singleM1), 1, fx)## Warning in wilcox.test.default(x = structure(c(0.017655648, 0.032979609, :
## cannot compute exact p-value with ties
## Warning in wilcox.test.default(x = structure(c(0, 0.000104037,
## 0.010320597, : cannot compute exact p-value with ties
## Warning in wilcox.test.default(x = structure(c(0, 0, 0.000878349,
## 0.003024803, : cannot compute exact p-value with ties
## Warning in wilcox.test.default(x = structure(c(0, 0.00031211,
## 0.000109794, : cannot compute exact p-value with ties
## Warning in wilcox.test.default(x = structure(c(0, 0, 0.000548968,
## 0.000950652, : cannot compute exact p-value with ties
## Warning in wilcox.test.default(x = structure(c(0.000132749, 0.00062422, :
## cannot compute exact p-value with ties
## Warning in wilcox.test.default(x = structure(c(0.000132749, 0.000104037, :
## cannot compute exact p-value with ties
## Warning in wilcox.test.default(x = structure(c(0, 0.000416146,
## 0.000329381, : cannot compute exact p-value with ties
## Warning in wilcox.test.default(x = structure(c(0, 0.00031211,
## 0.000439174, : cannot compute exact p-value with ties
## Warning in wilcox.test.default(x = structure(c(0.000132749, 0,
## 0.000109794, : cannot compute exact p-value with ties
## Warning in wilcox.test.default(x = structure(c(0, 0, 0, 0.000432115,
## 0.00088476: cannot compute exact p-value with ties
## Warning in wilcox.test.default(x = structure(c(0, 0.000104037, 0,
## 0.000345692, : cannot compute exact p-value with ties
## Warning in wilcox.test.default(x = structure(c(0, 0, 0.000219587,
## 0.000345692, : cannot compute exact p-value with ties
## Warning in wilcox.test.default(x = structure(c(0, 0, 0.000658762,
## 0.000345692, : cannot compute exact p-value with ties
## Warning in wilcox.test.default(x = structure(c(0, 0.000208073,
## 0.000219587, : cannot compute exact p-value with ties
## Warning in wilcox.test.default(x = structure(c(0, 0.000104037,
## 0.000548968, : cannot compute exact p-value with ties
## Warning in wilcox.test.default(x = structure(c(0.000132749, 0, 0,
## 0.000345692, : cannot compute exact p-value with ties
## Warning in wilcox.test.default(x = structure(c(0, 0.000104037, 0,
## 0.000259269, : cannot compute exact p-value with ties
## Warning in wilcox.test.default(x = structure(c(0, 0.000416146,
## 0.000109794, : cannot compute exact p-value with ties
## Warning in wilcox.test.default(x = structure(c(0, 0, 0.000878349,
## 0.000259269, : cannot compute exact p-value with ties
## Warning in wilcox.test.default(x = structure(c(0, 0.000104037,
## 0.000439174, : cannot compute exact p-value with ties
## Warning in wilcox.test.default(x = structure(c(0.000398248, 0.000208073, :
## cannot compute exact p-value with ties
## Warning in wilcox.test.default(x = structure(c(0.000663746, 0.000104037, :
## cannot compute exact p-value with ties
## Warning in wilcox.test.default(x = structure(c(0, 0, 0.000878349,
## 0.000172846, : cannot compute exact p-value with ties
## Warning in wilcox.test.default(x = structure(c(0.000265498, 0,
## 0.000439174, : cannot compute exact p-value with ties
## Warning in wilcox.test.default(x = structure(c(0.000132749, 0, 0,
## 0.000172846, : cannot compute exact p-value with ties
## Warning in wilcox.test.default(x = structure(c(0, 0.00031211, 0,
## 8.64e-05, : cannot compute exact p-value with ties
## Warning in wilcox.test.default(x = structure(c(0.000530997, 0.000208073, :
## cannot compute exact p-value with ties
## Warning in wilcox.test.default(x = structure(c(0, 0.000520183,
## 0.000219587, : cannot compute exact p-value with ties
## Warning in wilcox.test.default(x = structure(c(0, 0, 0.000109794,
## 8.64e-05, : cannot compute exact p-value with ties
## Warning in wilcox.test.default(x = structure(c(0.000530997, 0,
## 0.000109794, : cannot compute exact p-value with ties
## Warning in wilcox.test.default(x = structure(c(0.000132749, 0.000104037, :
## cannot compute exact p-value with ties
## Warning in wilcox.test.default(x = structure(c(0, 0.00031211, 0,
## 8.64e-05, : cannot compute exact p-value with ties
## Warning in wilcox.test.default(x = structure(c(0, 0.000208073, 0,
## 8.64e-05, : cannot compute exact p-value with ties
## Warning in wilcox.test.default(x = structure(c(0.000265498, 0,
## 0.000768555, : cannot compute exact p-value with ties
## Warning in wilcox.test.default(x = structure(c(0.000663746, 0.000104037, :
## cannot compute exact p-value with ties
## Warning in wilcox.test.default(x = structure(c(0, 0, 0.000109794,
## 8.64e-05, : cannot compute exact p-value with ties
## Warning in wilcox.test.default(x = structure(c(0, 0.000104037,
## 0.000109794, : cannot compute exact p-value with ties
## Warning in wilcox.test.default(x = structure(c(0, 0, 0.000109794,
## 8.64e-05, : cannot compute exact p-value with ties
## Warning in wilcox.test.default(x = structure(c(0.000132749, 0,
## 0.000109794, : cannot compute exact p-value with ties
## Warning in wilcox.test.default(x = structure(c(0, 0.000104037, 0,
## 8.64e-05, : cannot compute exact p-value with ties
## Warning in wilcox.test.default(x = structure(c(0, 0.000104037, 0,
## 8.64e-05, : cannot compute exact p-value with ties
## Warning in wilcox.test.default(x = structure(c(0.000265498, 0.000104037, :
## cannot compute exact p-value with ties
## Warning in wilcox.test.default(x = structure(c(0.000132749, 0.000104037, :
## cannot compute exact p-value with ties
## Warning in wilcox.test.default(x = structure(c(0, 0, 0, 8.64e-05,
## 0.00022119: cannot compute exact p-value with ties
## Warning in wilcox.test.default(x = structure(c(0, 0.000104037, 0,
## 8.64e-05, : cannot compute exact p-value with ties
## Warning in wilcox.test.default(x = structure(c(0, 0, 0, 8.64e-05,
## 0.00022119: cannot compute exact p-value with ties
## Warning in wilcox.test.default(x = structure(c(0.000132749, 0.000520183, :
## cannot compute exact p-value with ties
## Warning in wilcox.test.default(x = structure(c(0.000265498, 0, 0,
## 8.64e-05, : cannot compute exact p-value with ties
## Warning in wilcox.test.default(x = structure(c(0, 0, 0, 0, 0), .Names =
## c("Lam_20", : cannot compute exact p-value with ties
## Warning in wilcox.test.default(x = structure(c(0.000265498, 0.000832293, :
## cannot compute exact p-value with ties
## Warning in wilcox.test.default(x = structure(c(0, 0, 0.000109794, 0,
## 0), .Names = c("Lam_20", : cannot compute exact p-value with ties
## Warning in wilcox.test.default(x = structure(c(0, 0, 0.000109794, 0,
## 0), .Names = c("Lam_20", : cannot compute exact p-value with ties
## Warning in wilcox.test.default(x = structure(c(0, 0, 0, 0,
## 0.00022119), .Names = c("Lam_20", : cannot compute exact p-value with ties
## Warning in wilcox.test.default(x = structure(c(0, 0, 0.000548968, 0,
## 0), .Names = c("Lam_20", : cannot compute exact p-value with ties
## Warning in wilcox.test.default(x = structure(c(0, 0, 0.000109794, 0,
## 0), .Names = c("Lam_20", : cannot compute exact p-value with ties
## Warning in wilcox.test.default(x = structure(c(0, 0, 0, 0, 0), .Names =
## c("Lam_20", : cannot compute exact p-value with ties
## Warning in wilcox.test.default(x = structure(c(0, 0.000104037, 0, 0,
## 0), .Names = c("Lam_20", : cannot compute exact p-value with ties
## Warning in wilcox.test.default(x = structure(c(0.000132749, 0.00031211, :
## cannot compute exact p-value with ties
## Warning in wilcox.test.default(x = structure(c(0, 0, 0, 0, 0), .Names =
## c("Lam_20", : cannot compute exact p-value with ties
## Warning in wilcox.test.default(x = structure(c(0, 0, 0, 0,
## 0.000110595), .Names = c("Lam_20", : cannot compute exact p-value with ties
## Warning in wilcox.test.default(x = structure(c(0, 0.000104037, 0, 0,
## 0), .Names = c("Lam_20", : cannot compute exact p-value with ties
## Warning in wilcox.test.default(x = structure(c(0.000132749, 0.000104037, :
## cannot compute exact p-value with ties
## Warning in wilcox.test.default(x = structure(c(0, 0.000832293, 0, 0,
## 0.000110595: cannot compute exact p-value with ties
## Warning in wilcox.test.default(x = structure(c(0.000132749, 0,
## 0.000329381, : cannot compute exact p-value with ties
## Warning in wilcox.test.default(x = structure(c(0, 0, 0.000219587, 0,
## 0.00022119: cannot compute exact p-value with ties
## Warning in wilcox.test.default(x = structure(c(0, 0, 0.000109794, 0,
## 0.000331785: cannot compute exact p-value with ties
## Warning in wilcox.test.default(x = structure(c(0.000265498, 0.000104037, :
## cannot compute exact p-value with ties
## Warning in wilcox.test.default(x = structure(c(0.000132749, 0, 0, 0,
## 0), .Names = c("Lam_20", : cannot compute exact p-value with ties
## Warning in wilcox.test.default(x = structure(c(0, 0, 0, 0,
## 0.000110595), .Names = c("Lam_20", : cannot compute exact p-value with ties
## Warning in wilcox.test.default(x = structure(c(0, 0, 0, 0, 0), .Names =
## c("Lam_20", : cannot compute exact p-value with ties
## Warning in wilcox.test.default(x = structure(c(0, 0.000104037, 0, 0,
## 0), .Names = c("Lam_20", : cannot compute exact p-value with ties
## Warning in wilcox.test.default(x = structure(c(0, 0, 0, 0, 0), .Names =
## c("Lam_20", : cannot compute exact p-value with ties
## Warning in wilcox.test.default(x = structure(c(0.000132749, 0, 0, 0,
## 0), .Names = c("Lam_20", : cannot compute exact p-value with ties
## Warning in wilcox.test.default(x = structure(c(0, 0, 0, 0, 0), .Names =
## c("Lam_20", : cannot compute exact p-value with ties
## Warning in wilcox.test.default(x = structure(c(0.000132749, 0, 0, 0,
## 0.000110595: cannot compute exact p-value with ties
## Warning in wilcox.test.default(x = structure(c(0, 0, 0, 0, 0), .Names =
## c("Lam_20", : cannot compute exact p-value with ties
## Warning in wilcox.test.default(x = structure(c(0, 0, 0, 0, 0), .Names =
## c("Lam_20", : cannot compute exact p-value with ties
## Warning in wilcox.test.default(x = structure(c(0, 0.000208073, 0, 0,
## 0), .Names = c("Lam_20", : cannot compute exact p-value with ties
## Warning in wilcox.test.default(x = structure(c(0.000132749, 0, 0, 0,
## 0.000110595: cannot compute exact p-value with ties
## Warning in wilcox.test.default(x = structure(c(0, 0, 0, 0, 0), .Names =
## c("Lam_20", : cannot compute exact p-value with ties
## Warning in wilcox.test.default(x = structure(c(0, 0, 0, 0, 0), .Names =
## c("Lam_20", : cannot compute exact p-value with ties
## Warning in wilcox.test.default(x = structure(c(0, 0, 0, 0, 0), .Names =
## c("Lam_20", : cannot compute exact p-value with ties
## Warning in wilcox.test.default(x = structure(c(0.000132749, 0.000104037, :
## cannot compute exact p-value with ties
genus_singleM.wilcox1 <- data.frame(genus_singleM.wilcox1 )
colnames(genus_singleM.wilcox1) <- c("p.value")
write.csv(genus_singleM.wilcox1,file="genus_wilcox0.05_Lam_T4.csv")
rownames(genus_singleM.wilcox1)[genus_singleM.wilcox1$p.value< 0.01]## [1] "Lactobacillus" "Turicibacter" NA NA
## [5] NA
#**,p<0.01, "Turicibacter" ; "Lactobacillus"
rownames(genus_singleM.wilcox1)[genus_singleM.wilcox1$p.value< 0.05]## [1] " CAG-1031" "Lactobacillus" "Turicibacter" "Eubacterium_F"
## [5] "UBA9475" "Neglecta" "Bacteroides" "Ruminococcus_A"
## [9] NA NA NA
#*,p<0.05, "Ruminococcus_A"; "Bacteroides" ;"Neglecta"; "UBA9475";"Eubacterium_F";" CAG-1031"
#Lam and T10
genus_singleM2 <- genus_singleM[,c(idx_Lam, idx_T10)]
clin.sub_genus_singleM2 <- data.frame(ID = 1:ncol(genus_singleM2), group = rep(c("Lam","T10"), each=5))
fx <- function(x) {
res <-wilcox.test(x ~clin.sub_genus_singleM2$group)
return(res$p.value)
}
genus_singleM.wilcox2 <-apply(as.matrix(genus_singleM2), 1, fx)## Warning in wilcox.test.default(x = structure(c(0, 0.000104037,
## 0.010320597, : cannot compute exact p-value with ties
## Warning in wilcox.test.default(x = structure(c(0, 0, 0.000878349,
## 0.003024803, : cannot compute exact p-value with ties
## Warning in wilcox.test.default(x = structure(c(0, 0.00031211,
## 0.000109794, : cannot compute exact p-value with ties
## Warning in wilcox.test.default(x = structure(c(0, 0, 0.000548968,
## 0.000950652, : cannot compute exact p-value with ties
## Warning in wilcox.test.default(x = structure(c(0.000132749, 0.00062422, :
## cannot compute exact p-value with ties
## Warning in wilcox.test.default(x = structure(c(0.000132749, 0.000104037, :
## cannot compute exact p-value with ties
## Warning in wilcox.test.default(x = structure(c(0.001061994, 0.000208073, :
## cannot compute exact p-value with ties
## Warning in wilcox.test.default(x = structure(c(0, 0.00031211,
## 0.000439174, : cannot compute exact p-value with ties
## Warning in wilcox.test.default(x = structure(c(0.000132749, 0,
## 0.000109794, : cannot compute exact p-value with ties
## Warning in wilcox.test.default(x = structure(c(0, 0, 0, 0.000432115,
## 0.00088476: cannot compute exact p-value with ties
## Warning in wilcox.test.default(x = structure(c(0, 0.000104037, 0,
## 0.000345692, : cannot compute exact p-value with ties
## Warning in wilcox.test.default(x = structure(c(0, 0, 0.000219587,
## 0.000345692, : cannot compute exact p-value with ties
## Warning in wilcox.test.default(x = structure(c(0, 0, 0.000658762,
## 0.000345692, : cannot compute exact p-value with ties
## Warning in wilcox.test.default(x = structure(c(0, 0.000208073,
## 0.000219587, : cannot compute exact p-value with ties
## Warning in wilcox.test.default(x = structure(c(0, 0.000104037,
## 0.000548968, : cannot compute exact p-value with ties
## Warning in wilcox.test.default(x = structure(c(0.000132749, 0, 0,
## 0.000345692, : cannot compute exact p-value with ties
## Warning in wilcox.test.default(x = structure(c(0, 0.000104037, 0,
## 0.000259269, : cannot compute exact p-value with ties
## Warning in wilcox.test.default(x = structure(c(0, 0.000416146,
## 0.000109794, : cannot compute exact p-value with ties
## Warning in wilcox.test.default(x = structure(c(0, 0, 0.000878349,
## 0.000259269, : cannot compute exact p-value with ties
## Warning in wilcox.test.default(x = structure(c(0.000663746, 0.000104037, :
## cannot compute exact p-value with ties
## Warning in wilcox.test.default(x = structure(c(0, 0, 0.000878349,
## 0.000172846, : cannot compute exact p-value with ties
## Warning in wilcox.test.default(x = structure(c(0.000265498, 0,
## 0.000439174, : cannot compute exact p-value with ties
## Warning in wilcox.test.default(x = structure(c(0, 0.000208073,
## 0.000329381, : cannot compute exact p-value with ties
## Warning in wilcox.test.default(x = structure(c(0.000132749, 0, 0,
## 0.000172846, : cannot compute exact p-value with ties
## Warning in wilcox.test.default(x = structure(c(0, 0.00031211, 0,
## 8.64e-05, : cannot compute exact p-value with ties
## Warning in wilcox.test.default(x = structure(c(0.000530997, 0.000208073, :
## cannot compute exact p-value with ties
## Warning in wilcox.test.default(x = structure(c(0, 0.000520183,
## 0.000219587, : cannot compute exact p-value with ties
## Warning in wilcox.test.default(x = structure(c(0, 0, 0.000109794,
## 8.64e-05, : cannot compute exact p-value with ties
## Warning in wilcox.test.default(x = structure(c(0.000530997, 0,
## 0.000109794, : cannot compute exact p-value with ties
## Warning in wilcox.test.default(x = structure(c(0.000132749, 0.000104037, :
## cannot compute exact p-value with ties
## Warning in wilcox.test.default(x = structure(c(0, 0.00031211, 0,
## 8.64e-05, : cannot compute exact p-value with ties
## Warning in wilcox.test.default(x = structure(c(0, 0.000208073, 0,
## 8.64e-05, : cannot compute exact p-value with ties
## Warning in wilcox.test.default(x = structure(c(0.000265498, 0,
## 0.000768555, : cannot compute exact p-value with ties
## Warning in wilcox.test.default(x = structure(c(0.000663746, 0.000104037, :
## cannot compute exact p-value with ties
## Warning in wilcox.test.default(x = structure(c(0, 0, 0.000109794,
## 8.64e-05, : cannot compute exact p-value with ties
## Warning in wilcox.test.default(x = structure(c(0, 0.000104037,
## 0.000109794, : cannot compute exact p-value with ties
## Warning in wilcox.test.default(x = structure(c(0, 0, 0.000109794,
## 8.64e-05, : cannot compute exact p-value with ties
## Warning in wilcox.test.default(x = structure(c(0.000132749, 0,
## 0.000109794, : cannot compute exact p-value with ties
## Warning in wilcox.test.default(x = structure(c(0, 0.000104037, 0,
## 8.64e-05, : cannot compute exact p-value with ties
## Warning in wilcox.test.default(x = structure(c(0, 0.000104037, 0,
## 8.64e-05, : cannot compute exact p-value with ties
## Warning in wilcox.test.default(x = structure(c(0.000265498, 0.000104037, :
## cannot compute exact p-value with ties
## Warning in wilcox.test.default(x = structure(c(0.000132749, 0.000104037, :
## cannot compute exact p-value with ties
## Warning in wilcox.test.default(x = structure(c(0, 0, 0, 8.64e-05,
## 0.00022119: cannot compute exact p-value with ties
## Warning in wilcox.test.default(x = structure(c(0, 0.000104037, 0,
## 8.64e-05, : cannot compute exact p-value with ties
## Warning in wilcox.test.default(x = structure(c(0, 0, 0, 8.64e-05,
## 0.00022119: cannot compute exact p-value with ties
## Warning in wilcox.test.default(x = structure(c(0.000132749, 0.000520183, :
## cannot compute exact p-value with ties
## Warning in wilcox.test.default(x = structure(c(0.000265498, 0, 0,
## 8.64e-05, : cannot compute exact p-value with ties
## Warning in wilcox.test.default(x = structure(c(0, 0, 0, 0, 0), .Names =
## c("Lam_20", : cannot compute exact p-value with ties
## Warning in wilcox.test.default(x = structure(c(0, 0, 0.000109794, 0,
## 0), .Names = c("Lam_20", : cannot compute exact p-value with ties
## Warning in wilcox.test.default(x = structure(c(0, 0, 0.000109794, 0,
## 0), .Names = c("Lam_20", : cannot compute exact p-value with ties
## Warning in wilcox.test.default(x = structure(c(0, 0, 0, 0,
## 0.00022119), .Names = c("Lam_20", : cannot compute exact p-value with ties
## Warning in wilcox.test.default(x = structure(c(0, 0, 0.000548968, 0,
## 0), .Names = c("Lam_20", : cannot compute exact p-value with ties
## Warning in wilcox.test.default(x = structure(c(0, 0, 0.000109794, 0,
## 0), .Names = c("Lam_20", : cannot compute exact p-value with ties
## Warning in wilcox.test.default(x = structure(c(0, 0, 0, 0, 0), .Names =
## c("Lam_20", : cannot compute exact p-value with ties
## Warning in wilcox.test.default(x = structure(c(0, 0.000104037, 0, 0,
## 0), .Names = c("Lam_20", : cannot compute exact p-value with ties
## Warning in wilcox.test.default(x = structure(c(0.000132749, 0.00031211, :
## cannot compute exact p-value with ties
## Warning in wilcox.test.default(x = structure(c(0, 0, 0, 0, 0), .Names =
## c("Lam_20", : cannot compute exact p-value with ties
## Warning in wilcox.test.default(x = structure(c(0, 0, 0, 0,
## 0.000110595), .Names = c("Lam_20", : cannot compute exact p-value with ties
## Warning in wilcox.test.default(x = structure(c(0, 0.000104037, 0, 0,
## 0), .Names = c("Lam_20", : cannot compute exact p-value with ties
## Warning in wilcox.test.default(x = structure(c(0.000132749, 0.000104037, :
## cannot compute exact p-value with ties
## Warning in wilcox.test.default(x = structure(c(0, 0.000832293, 0, 0,
## 0.000110595: cannot compute exact p-value with ties
## Warning in wilcox.test.default(x = structure(c(0.000132749, 0,
## 0.000329381, : cannot compute exact p-value with ties
## Warning in wilcox.test.default(x = structure(c(0, 0, 0.000219587, 0,
## 0.00022119: cannot compute exact p-value with ties
## Warning in wilcox.test.default(x = structure(c(0, 0, 0.000109794, 0,
## 0.000331785: cannot compute exact p-value with ties
## Warning in wilcox.test.default(x = structure(c(0.000265498, 0.000104037, :
## cannot compute exact p-value with ties
## Warning in wilcox.test.default(x = structure(c(0.000132749, 0, 0, 0,
## 0), .Names = c("Lam_20", : cannot compute exact p-value with ties
## Warning in wilcox.test.default(x = structure(c(0, 0, 0, 0,
## 0.000110595), .Names = c("Lam_20", : cannot compute exact p-value with ties
## Warning in wilcox.test.default(x = structure(c(0, 0, 0, 0, 0), .Names =
## c("Lam_20", : cannot compute exact p-value with ties
## Warning in wilcox.test.default(x = structure(c(0, 0.000104037, 0, 0,
## 0), .Names = c("Lam_20", : cannot compute exact p-value with ties
## Warning in wilcox.test.default(x = structure(c(0, 0, 0, 0, 0), .Names =
## c("Lam_20", : cannot compute exact p-value with ties
## Warning in wilcox.test.default(x = structure(c(0.000132749, 0, 0, 0,
## 0), .Names = c("Lam_20", : cannot compute exact p-value with ties
## Warning in wilcox.test.default(x = structure(c(0, 0, 0, 0, 0), .Names =
## c("Lam_20", : cannot compute exact p-value with ties
## Warning in wilcox.test.default(x = structure(c(0.000132749, 0, 0, 0,
## 0.000110595: cannot compute exact p-value with ties
## Warning in wilcox.test.default(x = structure(c(0, 0, 0, 0, 0), .Names =
## c("Lam_20", : cannot compute exact p-value with ties
## Warning in wilcox.test.default(x = structure(c(0, 0, 0, 0, 0), .Names =
## c("Lam_20", : cannot compute exact p-value with ties
## Warning in wilcox.test.default(x = structure(c(0, 0.000208073, 0, 0,
## 0), .Names = c("Lam_20", : cannot compute exact p-value with ties
## Warning in wilcox.test.default(x = structure(c(0.000132749, 0, 0, 0,
## 0.000110595: cannot compute exact p-value with ties
## Warning in wilcox.test.default(x = structure(c(0, 0, 0, 0, 0), .Names =
## c("Lam_20", : cannot compute exact p-value with ties
## Warning in wilcox.test.default(x = structure(c(0, 0, 0, 0, 0), .Names =
## c("Lam_20", : cannot compute exact p-value with ties
## Warning in wilcox.test.default(x = structure(c(0, 0, 0, 0, 0), .Names =
## c("Lam_20", : cannot compute exact p-value with ties
## Warning in wilcox.test.default(x = structure(c(0.000132749, 0.000104037, :
## cannot compute exact p-value with ties
genus_singleM.wilcox2 <- data.frame(genus_singleM.wilcox2)
colnames(genus_singleM.wilcox2) <- c("p.value")
write.csv(genus_singleM.wilcox2,file="genus_wilcox0.05_Lam_T10.csv")
rownames(genus_singleM.wilcox2)[genus_singleM.wilcox2$p.value< 0.01]## [1] "Eubacterium_R" NA NA
#**, p<0.01, "Eubacterium_R"
rownames(genus_singleM.wilcox2)[genus_singleM.wilcox2$p.value< 0.05]## [1] "Lactobacillus" "Turicibacter" "Eubacterium_R"
## [4] "Clostridium_M" "Lachnospira" "Flavonifractor"
## [7] "Ruminococcus" "Slackia" "Pseudooceanicola"
## [10] NA NA
#*, p<0.05, "Slackia";"Pseudooceanicola";"Ruminococcus";"Flavonifractor";"Lachnospira";"Clostridium_M";"Turicibacter";"Lactobacillus"
####T4 and T10
genus_singleM3 <- genus_singleM[,c(idx_T4, idx_T10)]
clin.sub_genus_singleM3 <- data.frame(ID = 1:ncol(genus_singleM3), group = rep(c("T4","T10"), each=5))
fx <- function(x) {
res <-wilcox.test(x ~clin.sub_genus_singleM3$group)
return(res$p.value)
}
genus_singleM.wilcox3 <-apply(as.matrix(genus_singleM3), 1, fx)## Warning in wilcox.test.default(x = structure(c(0.0081413, 0.00039267,
## 0.000274085, : cannot compute exact p-value with ties
## Warning in wilcox.test.default(x = structure(c(0.021112184, 0.021858639, :
## cannot compute exact p-value with ties
## Warning in wilcox.test.default(x = structure(c(0, 0, 0, 0, 0), .Names =
## c("T10_14", : cannot compute exact p-value with ties
## Warning in wilcox.test.default(x = structure(c(0, 0.00078534,
## 0.000685213, : cannot compute exact p-value with ties
## Warning in wilcox.test.default(x = structure(c(0.002483786, 0,
## 0.001507469, : cannot compute exact p-value with ties
## Warning in wilcox.test.default(x = structure(c(0.000413964, 0.00091623, :
## cannot compute exact p-value with ties
## Warning in wilcox.test.default(x = structure(c(0.000137988, 0.00104712, :
## cannot compute exact p-value with ties
## Warning in wilcox.test.default(x = structure(c(0, 0, 0, 0.000104899,
## 8.11e-05: cannot compute exact p-value with ties
## Warning in wilcox.test.default(x = structure(c(0.000137988, 0, 0, 0,
## 0), .Names = c("T10_14", : cannot compute exact p-value with ties
## Warning in wilcox.test.default(x = structure(c(0.000413964, 0.00026178, :
## cannot compute exact p-value with ties
## Warning in wilcox.test.default(x = structure(c(0, 0, 0, 0, 0), .Names =
## c("T10_14", : cannot compute exact p-value with ties
## Warning in wilcox.test.default(x = structure(c(0.000137988, 0,
## 0.000411128, : cannot compute exact p-value with ties
## Warning in wilcox.test.default(x = structure(c(0, 0, 0.000274085,
## 0.000209798, : cannot compute exact p-value with ties
## Warning in wilcox.test.default(x = structure(c(0.000689941, 0.00026178, :
## cannot compute exact p-value with ties
## Warning in wilcox.test.default(x = structure(c(0, 0.00039267, 0, 0,
## 0.000324254: cannot compute exact p-value with ties
## Warning in wilcox.test.default(x = structure(c(0, 0, 0, 0.000104899,
## 8.11e-05: cannot compute exact p-value with ties
## Warning in wilcox.test.default(x = structure(c(0, 0.00026178, 0.00054817, :
## cannot compute exact p-value with ties
## Warning in wilcox.test.default(x = structure(c(0.000137988, 0,
## 0.000411128, : cannot compute exact p-value with ties
## Warning in wilcox.test.default(x = structure(c(0, 0, 0.000274085, 0,
## 0), .Names = c("T10_14", : cannot compute exact p-value with ties
## Warning in wilcox.test.default(x = structure(c(0.000137988, 0.00013089, :
## cannot compute exact p-value with ties
## Warning in wilcox.test.default(x = structure(c(0.000275976, 0.00104712, :
## cannot compute exact p-value with ties
## Warning in wilcox.test.default(x = structure(c(0, 0, 0.000274085,
## 0.000209798, : cannot compute exact p-value with ties
## Warning in wilcox.test.default(x = structure(c(0, 0, 0, 0.000524494,
## 0), .Names = c("T10_14", : cannot compute exact p-value with ties
## Warning in wilcox.test.default(x = structure(c(0.001379881, 0.00091623, :
## cannot compute exact p-value with ties
## Warning in wilcox.test.default(x = structure(c(0, 0.00026178, 0,
## 0.000419595, : cannot compute exact p-value with ties
## Warning in wilcox.test.default(x = structure(c(0.000137988, 0, 0,
## 0.000314696, : cannot compute exact p-value with ties
## Warning in wilcox.test.default(x = structure(c(0, 0, 0, 0.000104899,
## 0), .Names = c("T10_14", : cannot compute exact p-value with ties
## Warning in wilcox.test.default(x = structure(c(0.000965917, 0.00052356, :
## cannot compute exact p-value with ties
## Warning in wilcox.test.default(x = structure(c(0, 0.00026178,
## 0.000685213, : cannot compute exact p-value with ties
## Warning in wilcox.test.default(x = structure(c(0.000689941, 0.001570681, :
## cannot compute exact p-value with ties
## Warning in wilcox.test.default(x = structure(c(0, 0, 0.000137043, 0,
## 0), .Names = c("T10_14", : cannot compute exact p-value with ties
## Warning in wilcox.test.default(x = structure(c(0.000689941, 0,
## 0.000137043, : cannot compute exact p-value with ties
## Warning in wilcox.test.default(x = structure(c(0, 0, 0.000137043,
## 0.000104899, : cannot compute exact p-value with ties
## Warning in wilcox.test.default(x = structure(c(0, 0, 0.000137043, 0,
## 8.11e-05: cannot compute exact p-value with ties
## Warning in wilcox.test.default(x = structure(c(0.000137988, 0,
## 0.000137043, : cannot compute exact p-value with ties
## Warning in wilcox.test.default(x = structure(c(0.000137988, 0.001439791, :
## cannot compute exact p-value with ties
## Warning in wilcox.test.default(x = structure(c(0.000827929, 0.00013089, :
## cannot compute exact p-value with ties
## Warning in wilcox.test.default(x = structure(c(0.000137988, 0, 0, 0,
## 0), .Names = c("T10_14", : cannot compute exact p-value with ties
## Warning in wilcox.test.default(x = structure(c(0, 0, 0, 0, 0), .Names =
## c("T10_14", : cannot compute exact p-value with ties
## Warning in wilcox.test.default(x = structure(c(0, 0, 0, 0, 0), .Names =
## c("T10_14", : cannot compute exact p-value with ties
## Warning in wilcox.test.default(x = structure(c(0.000413964, 0.00013089, :
## cannot compute exact p-value with ties
## Warning in wilcox.test.default(x = structure(c(0, 0, 0, 0.000104899,
## 0), .Names = c("T10_14", : cannot compute exact p-value with ties
## Warning in wilcox.test.default(x = structure(c(0, 0.00052356, 0,
## 0.000104899, : cannot compute exact p-value with ties
## Warning in wilcox.test.default(x = structure(c(0, 0.00013089, 0, 0,
## 8.11e-05: cannot compute exact p-value with ties
## Warning in wilcox.test.default(x = structure(c(0, 0, 0, 0.000104899,
## 8.11e-05: cannot compute exact p-value with ties
## Warning in wilcox.test.default(x = structure(c(0, 0.00013089, 0, 0,
## 8.11e-05: cannot compute exact p-value with ties
## Warning in wilcox.test.default(x = structure(c(0, 0, 0, 0, 0), .Names =
## c("T10_14", : cannot compute exact p-value with ties
## Warning in wilcox.test.default(x = structure(c(0, 0, 0, 0, 0), .Names =
## c("T10_14", : cannot compute exact p-value with ties
## Warning in wilcox.test.default(x = structure(c(0, 0.00013089, 0,
## 0.000524494, : cannot compute exact p-value with ties
## Warning in wilcox.test.default(x = structure(c(0, 0, 0, 0, 0), .Names =
## c("T10_14", : cannot compute exact p-value with ties
## Warning in wilcox.test.default(x = structure(c(0, 0.00013089,
## 0.001644511, : cannot compute exact p-value with ties
## Warning in wilcox.test.default(x = structure(c(0, 0.00026178,
## 0.000274085, : cannot compute exact p-value with ties
## Warning in wilcox.test.default(x = structure(c(0, 0.00013089,
## 0.000274085, : cannot compute exact p-value with ties
## Warning in wilcox.test.default(x = structure(c(0, 0.00026178,
## 0.000274085, : cannot compute exact p-value with ties
## Warning in wilcox.test.default(x = structure(c(0.000137988, 0.00078534, :
## cannot compute exact p-value with ties
## Warning in wilcox.test.default(x = structure(c(0, 0.00013089,
## 0.000137043, : cannot compute exact p-value with ties
## Warning in wilcox.test.default(x = structure(c(0, 0.00013089,
## 0.000137043, : cannot compute exact p-value with ties
## Warning in wilcox.test.default(x = structure(c(0.000137988, 0,
## 0.000137043, : cannot compute exact p-value with ties
## Warning in wilcox.test.default(x = structure(c(0, 0, 0.000137043,
## 0.000104899, : cannot compute exact p-value with ties
## Warning in wilcox.test.default(x = structure(c(0, 0.00026178,
## 0.000137043, : cannot compute exact p-value with ties
## Warning in wilcox.test.default(x = structure(c(0, 0.00013089,
## 0.000137043, : cannot compute exact p-value with ties
## Warning in wilcox.test.default(x = structure(c(0.000137988, 0,
## 0.000137043, : cannot compute exact p-value with ties
## Warning in wilcox.test.default(x = structure(c(0, 0.00013089,
## 0.000137043, : cannot compute exact p-value with ties
## Warning in wilcox.test.default(x = structure(c(0, 0, 0.000137043, 0,
## 0), .Names = c("T10_14", : cannot compute exact p-value with ties
## Warning in wilcox.test.default(x = structure(c(0, 0, 0, 0,
## 8.11e-05), .Names = c("T10_14", : cannot compute exact p-value with ties
## Warning in wilcox.test.default(x = structure(c(0.000551953, 0.00013089, :
## cannot compute exact p-value with ties
## Warning in wilcox.test.default(x = structure(c(0.000137988, 0, 0, 0,
## 0), .Names = c("T10_14", : cannot compute exact p-value with ties
## Warning in wilcox.test.default(x = structure(c(0, 0, 0, 0,
## 8.11e-05), .Names = c("T10_14", : cannot compute exact p-value with ties
## Warning in wilcox.test.default(x = structure(c(0, 0.00013089, 0, 0,
## 0), .Names = c("T10_14", : cannot compute exact p-value with ties
## Warning in wilcox.test.default(x = structure(c(0.000137988, 0, 0, 0,
## 8.11e-05: cannot compute exact p-value with ties
## Warning in wilcox.test.default(x = structure(c(0, 0, 0, 0, 0), .Names =
## c("T10_14", : cannot compute exact p-value with ties
## Warning in wilcox.test.default(x = structure(c(0, 0, 0, 0.000104899,
## 8.11e-05: cannot compute exact p-value with ties
## Warning in wilcox.test.default(x = structure(c(0, 0.00026178, 0, 0,
## 0), .Names = c("T10_14", : cannot compute exact p-value with ties
## Warning in wilcox.test.default(x = structure(c(0.000965917, 0, 0,
## 0.000209798, : cannot compute exact p-value with ties
## Warning in wilcox.test.default(x = structure(c(0.000137988, 0, 0,
## 0.000104899, : cannot compute exact p-value with ties
## Warning in wilcox.test.default(x = structure(c(0, 0, 0, 0, 0), .Names =
## c("T10_14", : cannot compute exact p-value with ties
## Warning in wilcox.test.default(x = structure(c(0, 0, 0, 0, 0), .Names =
## c("T10_14", : cannot compute exact p-value with ties
## Warning in wilcox.test.default(x = structure(c(0.000137988, 0.00013089, :
## cannot compute exact p-value with ties
## Warning in wilcox.test.default(x = structure(c(0, 0.00026178, 0, 0,
## 8.11e-05: cannot compute exact p-value with ties
## Warning in wilcox.test.default(x = structure(c(0, 0, 0, 0, 0), .Names =
## c("T10_14", : cannot compute exact p-value with ties
## Warning in wilcox.test.default(x = structure(c(0, 0, 0, 0.000209798,
## 0), .Names = c("T10_14", : cannot compute exact p-value with ties
## Warning in wilcox.test.default(x = structure(c(0, 0, 0, 0.000104899,
## 0), .Names = c("T10_14", : cannot compute exact p-value with ties
## Warning in wilcox.test.default(x = structure(c(0, 0.00013089, 0,
## 0.000209798, : cannot compute exact p-value with ties
## Warning in wilcox.test.default(x = structure(c(0, 0, 0, 0,
## 8.11e-05), .Names = c("T10_14", : cannot compute exact p-value with ties
genus_singleM.wilcox3 <- data.frame(genus_singleM.wilcox3)
colnames(genus_singleM.wilcox3) <- c("p.value")
write.csv(genus_singleM.wilcox3,file="genus_wilcox0.05_T4_T10.csv")
rownames(genus_singleM.wilcox3)[genus_singleM.wilcox3$p.value< 0.01]## [1] "Eubacterium_I" NA NA NA
## [5] "CAG-791" NA
#**,p<0.01, "CAG-791","Eubacterium_I"
rownames(genus_singleM.wilcox3)[genus_singleM.wilcox3$p.value< 0.05]## [1] "Clostridium_M" "Eubacterium_I" "Weissella" NA
## [5] NA NA "CAG-791" NA
#*,p<0.05, "Weissella","Clostridium_M"#pairwise.wilcox.test with fdr correction
genus_singleM4=genus_singleM*100
genus_singleM5=t(genus_singleM4)
singleM_genus=cbind(genus_singleM5,Colvec_expedition)
pairwise.wilcox.test(singleM_genus$` CAG-1031`, singleM_genus$Zones,p.adjust.method = "fdr")##
## Pairwise comparisons using Wilcoxon rank sum test
##
## data: singleM_genus$` CAG-1031` and singleM_genus$Zones
##
## Lam T10
## T10 0.143 -
## T4 0.048 0.690
##
## P value adjustment method: fdr
pairwise.wilcox.test(singleM_genus$Lactobacillus, singleM_genus$Zones,p.adjust.method = "fdr")##
## Pairwise comparisons using Wilcoxon rank sum test
##
## data: singleM_genus$Lactobacillus and singleM_genus$Zones
##
## Lam T10
## T10 0.048 -
## T4 0.024 0.841
##
## P value adjustment method: fdr
pairwise.wilcox.test(singleM_genus$Turicibacter, singleM_genus$Zones,p.adjust.method = "fdr")## Warning in wilcox.test.default(xi, xj, paired = paired, ...): cannot
## compute exact p-value with ties
## Warning in wilcox.test.default(xi, xj, paired = paired, ...): cannot
## compute exact p-value with ties
##
## Pairwise comparisons using Wilcoxon rank sum test
##
## data: singleM_genus$Turicibacter and singleM_genus$Zones
##
## Lam T10
## T10 0.048 -
## T4 0.029 0.181
##
## P value adjustment method: fdr
pairwise.wilcox.test(singleM_genus$Bacteroides, singleM_genus$Zones,p.adjust.method = "fdr")##
## Pairwise comparisons using Wilcoxon rank sum test
##
## data: singleM_genus$Bacteroides and singleM_genus$Zones
##
## Lam T10
## T10 0.690 -
## T4 0.048 0.690
##
## P value adjustment method: fdr
pairwise.wilcox.test(singleM_genus$Weissella, singleM_genus$Zones,p.adjust.method = "fdr")## Warning in wilcox.test.default(xi, xj, paired = paired, ...): cannot
## compute exact p-value with ties
## Warning in wilcox.test.default(xi, xj, paired = paired, ...): cannot
## compute exact p-value with ties
## Warning in wilcox.test.default(xi, xj, paired = paired, ...): cannot
## compute exact p-value with ties
##
## Pairwise comparisons using Wilcoxon rank sum test
##
## data: singleM_genus$Weissella and singleM_genus$Zones
##
## Lam T10
## T10 0.8 -
## T4 0.1 0.1
##
## P value adjustment method: fdr
#members in class Clostridia
pairwise.wilcox.test(singleM_genus$Eubacterium_R, singleM_genus$Zones,p.adjust.method = "fdr")##
## Pairwise comparisons using Wilcoxon rank sum test
##
## data: singleM_genus$Eubacterium_R and singleM_genus$Zones
##
## Lam T10
## T10 0.024 -
## T4 0.548 0.226
##
## P value adjustment method: fdr
pairwise.wilcox.test(singleM_genus$Eubacterium_F, singleM_genus$Zones,p.adjust.method = "fdr")##
## Pairwise comparisons using Wilcoxon rank sum test
##
## data: singleM_genus$Eubacterium_F and singleM_genus$Zones
##
## Lam T10
## T10 0.226 -
## T4 0.048 0.841
##
## P value adjustment method: fdr
pairwise.wilcox.test(singleM_genus$UBA9475, singleM_genus$Zones,p.adjust.method = "fdr")##
## Pairwise comparisons using Wilcoxon rank sum test
##
## data: singleM_genus$UBA9475 and singleM_genus$Zones
##
## Lam T10
## T10 0.690 -
## T4 0.095 0.143
##
## P value adjustment method: fdr
pairwise.wilcox.test(singleM_genus$Neglecta, singleM_genus$Zones,p.adjust.method = "fdr")##
## Pairwise comparisons using Wilcoxon rank sum test
##
## data: singleM_genus$Neglecta and singleM_genus$Zones
##
## Lam T10
## T10 0.548 -
## T4 0.095 0.464
##
## P value adjustment method: fdr
pairwise.wilcox.test(singleM_genus$Lachnospira, singleM_genus$Zones,p.adjust.method = "fdr")##
## Pairwise comparisons using Wilcoxon rank sum test
##
## data: singleM_genus$Lachnospira and singleM_genus$Zones
##
## Lam T10
## T10 0.095 -
## T4 0.143 0.421
##
## P value adjustment method: fdr
#rare genera
pairwise.wilcox.test(singleM_genus$Clostridium_M, singleM_genus$Zones,p.adjust.method = "fdr")##
## Pairwise comparisons using Wilcoxon rank sum test
##
## data: singleM_genus$Clostridium_M and singleM_genus$Zones
##
## Lam T10
## T10 0.024 -
## T4 1.000 0.024
##
## P value adjustment method: fdr
pairwise.wilcox.test(singleM_genus$Eubacterium_I, singleM_genus$Zones,p.adjust.method = "fdr")##
## Pairwise comparisons using Wilcoxon rank sum test
##
## data: singleM_genus$Eubacterium_I and singleM_genus$Zones
##
## Lam T10
## T10 0.222 -
## T4 0.222 0.024
##
## P value adjustment method: fdr
pairwise.wilcox.test(singleM_genus$`CAG-791`, singleM_genus$Zones,p.adjust.method = "fdr")## Warning in wilcox.test.default(xi, xj, paired = paired, ...): cannot
## compute exact p-value with ties
## Warning in wilcox.test.default(xi, xj, paired = paired, ...): cannot
## compute exact p-value with ties
## Warning in wilcox.test.default(xi, xj, paired = paired, ...): cannot
## compute exact p-value with ties
##
## Pairwise comparisons using Wilcoxon rank sum test
##
## data: singleM_genus$`CAG-791` and singleM_genus$Zones
##
## Lam T10
## T10 0.441 -
## T4 0.086 0.029
##
## P value adjustment method: fdr
pairwise.wilcox.test(singleM_genus$Slackia, singleM_genus$Zones,p.adjust.method = "fdr")## Warning in wilcox.test.default(xi, xj, paired = paired, ...): cannot
## compute exact p-value with ties
## Warning in wilcox.test.default(xi, xj, paired = paired, ...): cannot
## compute exact p-value with ties
## Warning in wilcox.test.default(xi, xj, paired = paired, ...): cannot
## compute exact p-value with ties
##
## Pairwise comparisons using Wilcoxon rank sum test
##
## data: singleM_genus$Slackia and singleM_genus$Zones
##
## Lam T10
## T10 0.033 -
## T4 0.373 0.141
##
## P value adjustment method: fdr
pairwise.wilcox.test(singleM_genus$Pseudooceanicola, singleM_genus$Zones,p.adjust.method = "fdr")## Warning in wilcox.test.default(xi, xj, paired = paired, ...): cannot
## compute exact p-value with ties
## Warning in wilcox.test.default(xi, xj, paired = paired, ...): cannot
## compute exact p-value with ties
## Warning in wilcox.test.default(xi, xj, paired = paired, ...): cannot
## compute exact p-value with ties
##
## Pairwise comparisons using Wilcoxon rank sum test
##
## data: singleM_genus$Pseudooceanicola and singleM_genus$Zones
##
## Lam T10
## T10 0.076 -
## T4 0.914 0.270
##
## P value adjustment method: fdr
pairwise.wilcox.test(singleM_genus$Ruminococcus, singleM_genus$Zones,p.adjust.method = "fdr")## Warning in wilcox.test.default(xi, xj, paired = paired, ...): cannot
## compute exact p-value with ties
## Warning in wilcox.test.default(xi, xj, paired = paired, ...): cannot
## compute exact p-value with ties
## Warning in wilcox.test.default(xi, xj, paired = paired, ...): cannot
## compute exact p-value with ties
##
## Pairwise comparisons using Wilcoxon rank sum test
##
## data: singleM_genus$Ruminococcus and singleM_genus$Zones
##
## Lam T10
## T10 0.076 -
## T4 0.833 0.108
##
## P value adjustment method: fdr
pairwise.wilcox.test(singleM_genus$Flavonifractor, singleM_genus$Zones,p.adjust.method = "fdr")##
## Pairwise comparisons using Wilcoxon rank sum test
##
## data: singleM_genus$Flavonifractor and singleM_genus$Zones
##
## Lam T10
## T10 0.048 -
## T4 0.151 0.143
##
## P value adjustment method: fdr
pairwise.wilcox.test(singleM_genus$Ruminococcus_A, singleM_genus$Zones,p.adjust.method = "fdr")## Warning in wilcox.test.default(xi, xj, paired = paired, ...): cannot
## compute exact p-value with ties
## Warning in wilcox.test.default(xi, xj, paired = paired, ...): cannot
## compute exact p-value with ties
##
## Pairwise comparisons using Wilcoxon rank sum test
##
## data: singleM_genus$Ruminococcus_A and singleM_genus$Zones
##
## Lam T10
## T10 0.42 -
## T4 0.11 0.42
##
## P value adjustment method: fdr
#Box-plots for significant genus
#A. CAG-1031
CAG=ggplot(singleM_genus, aes(x=Zones, y=` CAG-1031`)) + geom_boxplot(fill=c("#999999", "#E69F00", "#56B4E9"))+geom_dotplot(binaxis='y', stackdir='center', dotsize=0.6)+
labs(title = NULL,x=NULL,y = "CAG-1031 \n relative abundances(%)")+theme_classic()
CAG_boxplot=CAG+theme(plot.title = element_text(size = 28,hjust=0.5))+theme(axis.text.x = element_text(color="black",size = 20))+theme(axis.text.y = element_text(color="black",size = 18),axis.title.y=element_text(size=20))
CAG_boxplot+scale_y_continuous(expand = c(0,0), limits = c(0,40))+geom_segment(y=36,yend=36,x=1,xend=3)+geom_text(x=2,y=36,label="*",size=12)+theme(plot.margin=unit(c(1,1,1.5,1.2),"cm"))+theme(axis.title.y = element_text(margin = margin(t = 0, r = 20, b = 0, l = 0)))## `stat_bindot()` using `bins = 30`. Pick better value with `binwidth`.
#B. Lactobacillus
Lactobacillus=ggplot(singleM_genus, aes(x=Zones, y=Lactobacillus)) + geom_boxplot(fill=c("#999999", "#E69F00", "#56B4E9"))+geom_dotplot(binaxis='y', stackdir='center', dotsize=0.6)+
labs(title = NULL,x=NULL,y = "Lactobacillus \n relative abundances(%)")+theme_classic()
Lactobacillus_boxplot=Lactobacillus+theme(plot.title = element_text(size = 28,hjust=0.5))+theme(axis.text.x = element_text(color="black",size = 20))+theme(axis.text.y = element_text(color="black",size = 18),axis.title.y=element_text(size=20))
Lactobacillus_boxplot+scale_y_continuous(expand = c(0,0), limits = c(0,10))+geom_segment(y=8.8,yend=8.8,x=1,xend=2)+geom_segment(y=9.5,yend=9.5,x=1,xend=3)+geom_text(x=1.5,y=8.8,label="*",size=12)+geom_text(x=2,y=9.5,label="**",size=12)+theme(plot.margin=unit(c(1,1,1.5,1.2),"cm"))+theme(axis.title.y = element_text(margin = margin(t = 0, r = 20, b = 0, l = 0)))## `stat_bindot()` using `bins = 30`. Pick better value with `binwidth`.
#Turicibacter
Turicibacter=ggplot(singleM_genus, aes(x=Zones, y=Turicibacter)) + geom_boxplot(fill=c("#999999", "#E69F00", "#56B4E9"))+geom_dotplot(binaxis='y', stackdir='center', dotsize=0.6)+
labs(title = NULL,x=NULL,y = "Turicibacter \n relative abundances(%)")+theme_classic()
Turicibacter_boxplot=Turicibacter+theme(plot.title = element_text(size = 28,hjust=0.5))+theme(axis.text.x = element_text(color="black",size = 20))+theme(axis.text.y = element_text(color="black",size = 18),axis.title.y=element_text(size=20))
Turicibacter_boxplot+scale_y_continuous(expand = c(0,0), limits = c(0,4.6))+geom_segment(y=3.8,yend=3.8,x=1,xend=2)+geom_segment(y=4.2,yend=4.2,x=1,xend=3)+geom_text(x=1.5,y=3.8,label="*",size=12)+geom_text(x=2,y=4.2,label="**",size=12)+theme(plot.margin=unit(c(1,1,1.5,1.2),"cm"))+theme(axis.title.y = element_text(margin = margin(t = 0, r = 20, b = 0, l = 0)))## `stat_bindot()` using `bins = 30`. Pick better value with `binwidth`.
#D. Weissella
Weissella=ggplot(singleM_genus, aes(x=Zones, y=Weissella)) + geom_boxplot(fill=c("#999999", "#E69F00", "#56B4E9"))+geom_dotplot(binaxis='y', stackdir='center', dotsize=0.6)+
labs(title =NULL,x=NULL,y = "Weissella \n relative abundances(%)")+theme_classic()
Weissella_boxplot=Weissella+theme(plot.title = element_text(size = 28,hjust=0.5))+theme(axis.text.x = element_text(color="black",size = 20))+theme(axis.text.y = element_text(color="black",size = 18),axis.title.y=element_text(size=20))
Weissella_boxplot+scale_y_continuous(expand = c(0,0), limits = c(0,12.5))+geom_segment(y=10.5,yend=10.5,x=2,xend=3)+geom_segment(y=11.2,yend=11.2,x=1,xend=3)+geom_text(x=2,y=11.8,label="p=0.07",size=8)+geom_text(x=2.5,y=10.5,label="*",size=12)+theme(plot.margin=unit(c(1,1,1.5,1.2),"cm"))+theme(axis.title.y = element_text(margin = margin(t = 0, r = 20, b = 0, l = 0)))## `stat_bindot()` using `bins = 30`. Pick better value with `binwidth`.
#E. Bacteroides
Bacteroides=ggplot(singleM_genus, aes(x=Zones, y=Bacteroides)) + geom_boxplot(fill=c("#999999", "#E69F00", "#56B4E9"))+geom_dotplot(binaxis='y', stackdir='center', dotsize=0.6)+
labs(title = NULL,x=NULL,y = "Bacteroides \n relative abundances(%)")+theme_classic()
Bacteroides_boxplot=Bacteroides+theme(plot.title = element_text(size = 28,hjust=0.5))+theme(axis.text.x = element_text(color="black",size = 20))+theme(axis.text.y = element_text(color="black",size = 18),axis.title.y=element_text(size=20))
Bacteroides_boxplot+scale_y_continuous(expand = c(0,0), limits = c(0,7))+geom_segment(y=6.4,yend=6.4,x=1,xend=3)+geom_text(x=2,y=6.4,label="*",size=12)+theme(plot.margin=unit(c(1,1,1.5,1.2),"cm"))+theme(axis.title.y = element_text(margin = margin(t = 0, r = 20, b = 0, l = 0)))## `stat_bindot()` using `bins = 30`. Pick better value with `binwidth`.
singleM_genus_Clostridia=data_frame(singleM_genus$Eubacterium_R,singleM_genus$Eubacterium_F,singleM_genus$UBA9475,singleM_genus$Neglecta,singleM_genus$Lachnospira,singleM_genus$Zones)## Warning: `data_frame()` is deprecated, use `tibble()`.
## This warning is displayed once per session.
colnames(singleM_genus_Clostridia)=c("Eubacterium_R","Eubacterium_F","UBA9475","Neglecta","Lachnospira","Zones")
row.names(singleM_genus_Clostridia)=row.names(singleM_genus)## Warning: Setting row names on a tibble is deprecated.
singleM_genus_Clostridia1=gather(singleM_genus_Clostridia,key="Clostridia",value = "abundances",-Zones)
genus_Clostridia_boxplot=ggplot(singleM_genus_Clostridia1, aes(x=Clostridia, y=abundances,fill=factor(Zones)))+theme_classic()+scale_fill_manual(values=c("#999999", "#E69F00", "#56B4E9"))+geom_boxplot(position=position_dodge(0.8))+ geom_dotplot(binaxis='y', stackdir='center', position=position_dodge(0.8),binpositions = "bygroup", dotsize=0.4)+labs(title =NULL,x=NULL,y = "Class Clostridia \n relative abundances(%)")+theme_classic()
genus_Clostridia_boxplot=genus_Clostridia_boxplot+theme(plot.title = element_text(size = 28,hjust=0.5))+theme(axis.text.x = element_text(color="black",size = 10,angle=45,hjust = 1))+theme(axis.text.y = element_text(color="black",size = 10),axis.title.y=element_text(size=12))
genus_Clostridia_boxplot+scale_y_continuous(expand = c(0,0), limits = c(0,4.5))+geom_text(x=1.275,y=2.2,label="*",size=8)+geom_text(x=2,y=3.6,label="*",size=8)+geom_text(x=3,y=1.1,label="*",size=8)+geom_text(x=4.275,y=2.2,label="*",size=8)+geom_text(x=5.275,y=1.2,label="*",size=8)+theme(axis.title.y = element_text(margin = margin(t = 0, r = 20, b = 0, l = 0)))## `stat_bindot()` using `bins = 30`. Pick better value with `binwidth`.
#read MAGs
mOTUs=read.table("microbial_bins_60percent_complete.txt",header=T,row.names=1,sep="\t")
mOTUs$Maxbins=as.factor(rownames((mOTUs)))
mOTUs_species=read.table("MAGs_species.txt",header=T,row.names=1,sep="\t")
mOTUs_species$Maxbins=as.factor(rownames((mOTUs_species)))
mOTUs_species1=merge(mOTUs_species,mOTUs,by="Maxbins")
mOTUs_species2=mOTUs_species1[,-1]
row.names(mOTUs_species2)=mOTUs_species2[,1]
mOTUs_species3=mOTUs_species2[,-1]
#Lam and T4
MAGs_species1 <- mOTUs_species3[,c(idx_Lam, idx_T4)]
clin.sub_MAGs_species1 <- data.frame(ID = 1:ncol(MAGs_species1), group = rep(c("Lam","T4"), each=5))
fx <- function(x) {
res <-wilcox.test(x ~clin.sub_MAGs_species1$group)
return(res$p.value)
}
MAGs_species.wilcox1 <-apply(as.matrix(MAGs_species1), 1, fx)## Warning in wilcox.test.default(x = structure(c(0.001, 0, 0, 0,
## 0.001), .Names = c("Lam_20", : cannot compute exact p-value with ties
## Warning in wilcox.test.default(x = structure(c(0, 0, 0, 0, 0), .Names =
## c("Lam_20", : cannot compute exact p-value with ties
## Warning in wilcox.test.default(x = structure(c(0, 0.001, 0.001, 0.001, 0:
## cannot compute exact p-value with ties
## Warning in wilcox.test.default(x = structure(c(0.214, 0.238, 0.248,
## 0.794, : cannot compute exact p-value with ties
## Warning in wilcox.test.default(x = structure(c(0.171, 0.188, 0.258,
## 0.537, : cannot compute exact p-value with ties
## Warning in wilcox.test.default(x = structure(c(0.226, 0.091, 0.107,
## 0.208, : cannot compute exact p-value with ties
## Warning in wilcox.test.default(x = structure(c(0.419, 0.426, 0.222,
## 0.276, : cannot compute exact p-value with ties
## Warning in wilcox.test.default(x = structure(c(0.011, 0.01, 0.023, 0.018, :
## cannot compute exact p-value with ties
## Warning in wilcox.test.default(x = structure(c(1.695, 3.866, 0.624,
## 0.248, : cannot compute exact p-value with ties
## Warning in wilcox.test.default(x = structure(c(0.035, 0.06, 0.006, 0.016, :
## cannot compute exact p-value with ties
## Warning in wilcox.test.default(x = structure(c(0.005, 0.007, 0.019,
## 0.016, : cannot compute exact p-value with ties
## Warning in wilcox.test.default(x = structure(c(0.633, 0.006, 0.22, 0.266, :
## cannot compute exact p-value with ties
## Warning in wilcox.test.default(x = structure(c(0.001, 0.002, 0.215,
## 0.935, : cannot compute exact p-value with ties
## Warning in wilcox.test.default(x = structure(c(0, 0, 1.775, 1.483, 3.445:
## cannot compute exact p-value with ties
MAGs_species.wilcox1 <- data.frame(MAGs_species.wilcox1 )
colnames(MAGs_species.wilcox1) <- c("p.value")
write.csv(MAGs_species.wilcox1, file = "MAGs_species.wilcox0.05_Lam_T4.csv")
rownames(MAGs_species.wilcox1)[MAGs_species.wilcox1$p.value< 0.01]## [1] "Lactococcus lactis_A" "maxbin.009"
## [3] "maxbin.020" "maxbin.026"
## [5] "maxbin.034" "Lactobacillus johnsonii"
## [7] "maxbin.087" "maxbin.121"
## [9] "maxbin.152"
#**, "Lactococcus lactis_A" "Lactobacillus johnsonii"
rownames(MAGs_species.wilcox1)[MAGs_species.wilcox1$p.value< 0.05]## [1] "Weissella cibaria" "Lactococcus lactis_A"
## [3] "maxbin.004" "maxbin.007"
## [5] "maxbin.009" "maxbin.014"
## [7] "maxbin.017" "maxbin.020"
## [9] "maxbin.026" "maxbin.029"
## [11] "maxbin.034" "Lactobacillus johnsonii"
## [13] "maxbin.057" "maxbin.062"
## [15] "maxbin.067" "maxbin.071"
## [17] "maxbin.087" "maxbin.121"
## [19] "maxbin.123" "maxbin.124"
## [21] "maxbin.133" "maxbin.152"
#*,"Weissella cibaria"
#Lam and T10
MAGs_species2 <- mOTUs_species3[,c(idx_Lam, idx_T10)]
clin.sub_MAGs_species2 <- data.frame(ID = 1:ncol(MAGs_species2), group = rep(c("Lam","T10"), each=5))
fx <- function(x) {
res <-wilcox.test(x ~clin.sub_MAGs_species2$group)
return(res$p.value)
}
MAGs_species.wilcox2 <-apply(as.matrix(MAGs_species2), 1, fx)## Warning in wilcox.test.default(x = structure(c(0.001, 0, 0, 0,
## 0.001), .Names = c("Lam_20", : cannot compute exact p-value with ties
## Warning in wilcox.test.default(x = structure(c(0, 0, 0, 0, 0), .Names =
## c("Lam_20", : cannot compute exact p-value with ties
## Warning in wilcox.test.default(x = structure(c(0.105, 0.127, 0.012,
## 0.052, : cannot compute exact p-value with ties
## Warning in wilcox.test.default(x = structure(c(0.617, 0.629, 0.501,
## 0.624, : cannot compute exact p-value with ties
## Warning in wilcox.test.default(x = structure(c(0, 0.001, 0.001, 0.001, 0:
## cannot compute exact p-value with ties
## Warning in wilcox.test.default(x = structure(c(0.171, 0.188, 0.258,
## 0.537, : cannot compute exact p-value with ties
## Warning in wilcox.test.default(x = structure(c(0.226, 0.091, 0.107,
## 0.208, : cannot compute exact p-value with ties
## Warning in wilcox.test.default(x = structure(c(0.04, 0.02, 4.768, 0.812, :
## cannot compute exact p-value with ties
## Warning in wilcox.test.default(x = structure(c(0.419, 0.426, 0.222,
## 0.276, : cannot compute exact p-value with ties
## Warning in wilcox.test.default(x = structure(c(1.207, 1.078, 0.041,
## 0.466, : cannot compute exact p-value with ties
## Warning in wilcox.test.default(x = structure(c(2.796, 5.897, 0.03,
## 18.361, : cannot compute exact p-value with ties
## Warning in wilcox.test.default(x = structure(c(2.227, 6.985, 0.006,
## 2.934, : cannot compute exact p-value with ties
## Warning in wilcox.test.default(x = structure(c(0.005, 0.007, 0.019,
## 0.016, : cannot compute exact p-value with ties
## Warning in wilcox.test.default(x = structure(c(0.004, 0.008, 0.412,
## 2.213, : cannot compute exact p-value with ties
## Warning in wilcox.test.default(x = structure(c(0.357, 0.241, 0.282,
## 0.537, : cannot compute exact p-value with ties
## Warning in wilcox.test.default(x = structure(c(0.001, 0.002, 0.215,
## 0.935, : cannot compute exact p-value with ties
## Warning in wilcox.test.default(x = structure(c(0, 0, 1.775, 1.483, 3.445:
## cannot compute exact p-value with ties
MAGs_species.wilcox2 <- data.frame(MAGs_species.wilcox2 )
colnames(MAGs_species.wilcox2) <- c("p.value")
write.csv(MAGs_species.wilcox2, file = "MAGs_species.wilcox0.05_Lam_T10.csv")
rownames(MAGs_species.wilcox2)[MAGs_species.wilcox2$p.value< 0.01]## [1] "Weissella cibaria" NA "maxbin.043"
## [4] "maxbin.087" "maxbin.123" "maxbin.141"
## [7] "maxbin.151" "maxbin.152" "maxbin.157"
## [10] "maxbin.166"
#**,"Weissella cibaria"
rownames(MAGs_species.wilcox2)[MAGs_species.wilcox2$p.value< 0.05]## [1] "Weissella cibaria" NA
## [3] "maxbin.011" "maxbin.014"
## [5] "maxbin.020" "maxbin.033"
## [7] "maxbin.043" "maxbin.044"
## [9] "Lactobacillus johnsonii" "maxbin.056"
## [11] "maxbin.087" "maxbin.121"
## [13] "maxbin.123" "maxbin.141"
## [15] "maxbin.151" "maxbin.152"
## [17] "maxbin.157" "maxbin.166"
#*,"Lactobacillus johnsonii"
#T10 and T4
MAGs_species3 <- mOTUs_species3[,c(idx_T4, idx_T10)]
clin.sub_MAGs_species3 <- data.frame(ID = 1:ncol(MAGs_species3), group = rep(c("T4","T10"), each=5))
fx <- function(x) {
res <-wilcox.test(x ~clin.sub_MAGs_species3$group)
return(res$p.value)
}
MAGs_species.wilcox3 <-apply(as.matrix(MAGs_species3), 1, fx)## Warning in wilcox.test.default(x = structure(c(0.005, 0.003, 0.002,
## 0.002, : cannot compute exact p-value with ties
## Warning in wilcox.test.default(x = structure(c(0, 0, 0, 0, 0), .Names =
## c("T10_14", : cannot compute exact p-value with ties
## Warning in wilcox.test.default(x = structure(c(0.21, 0.001, 0.001, 0.998, :
## cannot compute exact p-value with ties
## Warning in wilcox.test.default(x = structure(c(0.001, 0.001, 0.001,
## 0.001, : cannot compute exact p-value with ties
## Warning in wilcox.test.default(x = structure(c(0.015, 0.007, 0.012,
## 0.008, : cannot compute exact p-value with ties
## Warning in wilcox.test.default(x = structure(c(0.739, 0.89, 0.057, 0.228, :
## cannot compute exact p-value with ties
## Warning in wilcox.test.default(x = structure(c(0.049, 0.054, 0.08, 0.041, :
## cannot compute exact p-value with ties
## Warning in wilcox.test.default(x = structure(c(0.034, 0.017, 0.017,
## 0.079, : cannot compute exact p-value with ties
## Warning in wilcox.test.default(x = structure(c(5.504, 1.563, 0.02, 2.836, :
## cannot compute exact p-value with ties
## Warning in wilcox.test.default(x = structure(c(0.025, 0.014, 0.003,
## 0.095, : cannot compute exact p-value with ties
## Warning in wilcox.test.default(x = structure(c(0.045, 0.015, 14.483,
## 0.032, : cannot compute exact p-value with ties
## Warning in wilcox.test.default(x = structure(c(0.01, 0.015, 0.011, 0.013, :
## cannot compute exact p-value with ties
## Warning in wilcox.test.default(x = structure(c(0.027, 0.013, 0.003,
## 0.062, : cannot compute exact p-value with ties
## Warning in wilcox.test.default(x = structure(c(0.05, 0.044, 0.021, 4.964, :
## cannot compute exact p-value with ties
## Warning in wilcox.test.default(x = structure(c(0.01, 0.005, 0.01, 0.006, :
## cannot compute exact p-value with ties
## Warning in wilcox.test.default(x = structure(c(0.003, 0.003, 0.003,
## 0.002, : cannot compute exact p-value with ties
## Warning in wilcox.test.default(x = structure(c(0.295, 0.858, 0.001,
## 0.792, : cannot compute exact p-value with ties
## Warning in wilcox.test.default(x = structure(c(0.002, 0.004, 0.001,
## 0.003, : cannot compute exact p-value with ties
## Warning in wilcox.test.default(x = structure(c(6.245, 5.564, 0.005,
## 6.745, : cannot compute exact p-value with ties
MAGs_species.wilcox3 <- data.frame(MAGs_species.wilcox3 )
colnames(MAGs_species.wilcox3) <- c("p.value")
write.csv(MAGs_species.wilcox3, file = "MAGs_species.wilcox0.05_T4_T10.csv")
rownames(MAGs_species.wilcox3)[MAGs_species.wilcox3$p.value< 0.05]## [1] "Weissella cibaria" "Lactococcus lactis_A" "maxbin.019"
## [4] "maxbin.043" "maxbin.067" "maxbin.123"
## [7] "maxbin.135"
#**,"Lactococcus lactis_A"
rownames(MAGs_species.wilcox3)[MAGs_species.wilcox3$p.value< 0.05]## [1] "Weissella cibaria" "Lactococcus lactis_A" "maxbin.019"
## [4] "maxbin.043" "maxbin.067" "maxbin.123"
## [7] "maxbin.135"
#*, "Weissella cibaria" mOTUs_species4=t(mOTUs_species3)
mOTUs_species5=cbind(mOTUs_species4, Colvec_expedition)
#Wilcoxon Rank Sum test with 'fdr'('BH') correction for an FDR adjusted p-value (or q-value)
pairwise.wilcox.test(mOTUs_species5$`Lactobacillus johnsonii`, mOTUs_species5$Zones,p.adjust.method = "fdr")##
## Pairwise comparisons using Wilcoxon rank sum test
##
## data: mOTUs_species5$`Lactobacillus johnsonii` and mOTUs_species5$Zones
##
## Lam T10
## T10 0.048 -
## T4 0.024 0.841
##
## P value adjustment method: fdr
pairwise.wilcox.test(mOTUs_species5$`maxbin.009`, mOTUs_species5$Zones,p.adjust.method = "fdr")##
## Pairwise comparisons using Wilcoxon rank sum test
##
## data: mOTUs_species5$maxbin.009 and mOTUs_species5$Zones
##
## Lam T10
## T10 0.143 -
## T4 0.024 0.548
##
## P value adjustment method: fdr
pairwise.wilcox.test(mOTUs_species5$`maxbin.121`, mOTUs_species5$Zones,p.adjust.method = "fdr")##
## Pairwise comparisons using Wilcoxon rank sum test
##
## data: mOTUs_species5$maxbin.121 and mOTUs_species5$Zones
##
## Lam T10
## T10 0.024 -
## T4 0.024 0.548
##
## P value adjustment method: fdr
pairwise.wilcox.test(mOTUs_species5$`Weissella cibaria`, mOTUs_species5$Zones,p.adjust.method = "fdr")## Warning in wilcox.test.default(xi, xj, paired = paired, ...): cannot
## compute exact p-value with ties
## Warning in wilcox.test.default(xi, xj, paired = paired, ...): cannot
## compute exact p-value with ties
## Warning in wilcox.test.default(xi, xj, paired = paired, ...): cannot
## compute exact p-value with ties
##
## Pairwise comparisons using Wilcoxon rank sum test
##
## data: mOTUs_species5$`Weissella cibaria` and mOTUs_species5$Zones
##
## Lam T10
## T10 0.011 -
## T4 0.011 0.011
##
## P value adjustment method: fdr
pairwise.wilcox.test(mOTUs_species5$`Lactococcus lactis_A`, mOTUs_species5$Zones,p.adjust.method = "fdr")## Warning in wilcox.test.default(xi, xj, paired = paired, ...): cannot
## compute exact p-value with ties
## Warning in wilcox.test.default(xi, xj, paired = paired, ...): cannot
## compute exact p-value with ties
## Warning in wilcox.test.default(xi, xj, paired = paired, ...): cannot
## compute exact p-value with ties
##
## Pairwise comparisons using Wilcoxon rank sum test
##
## data: mOTUs_species5$`Lactococcus lactis_A` and mOTUs_species5$Zones
##
## Lam T10
## T10 - -
## T4 0.0075 0.0075
##
## P value adjustment method: fdr
pairwise.wilcox.test(mOTUs_species5$`Bacteroides thetaiotaomicron`, mOTUs_species5$Zones,p.adjust.method = "fdr")##
## Pairwise comparisons using Wilcoxon rank sum test
##
## data: mOTUs_species5$`Bacteroides thetaiotaomicron` and mOTUs_species5$Zones
##
## Lam T10
## T10 0.23 -
## T4 0.17 1.00
##
## P value adjustment method: fdr
pairwise.wilcox.test(mOTUs_species5$`maxbin.004`, mOTUs_species5$Zones,p.adjust.method = "fdr")##
## Pairwise comparisons using Wilcoxon rank sum test
##
## data: mOTUs_species5$maxbin.004 and mOTUs_species5$Zones
##
## Lam T10
## T10 0.083 -
## T4 0.083 0.151
##
## P value adjustment method: fdr
pairwise.wilcox.test(mOTUs_species5$`maxbin.007`, mOTUs_species5$Zones,p.adjust.method = "fdr")##
## Pairwise comparisons using Wilcoxon rank sum test
##
## data: mOTUs_species5$maxbin.007 and mOTUs_species5$Zones
##
## Lam T10
## T10 0.421 -
## T4 0.048 0.143
##
## P value adjustment method: fdr
pairwise.wilcox.test(mOTUs_species5$`maxbin.011`, mOTUs_species5$Zones,p.adjust.method = "fdr")##
## Pairwise comparisons using Wilcoxon rank sum test
##
## data: mOTUs_species5$maxbin.011 and mOTUs_species5$Zones
##
## Lam T10
## T10 0.095 -
## T4 0.226 1.000
##
## P value adjustment method: fdr
pairwise.wilcox.test(mOTUs_species5$`maxbin.014`, mOTUs_species5$Zones,p.adjust.method = "fdr")##
## Pairwise comparisons using Wilcoxon rank sum test
##
## data: mOTUs_species5$maxbin.014 and mOTUs_species5$Zones
##
## Lam T10
## T10 0.048 -
## T4 0.048 0.222
##
## P value adjustment method: fdr
pairwise.wilcox.test(mOTUs_species5$`maxbin.017`, mOTUs_species5$Zones,p.adjust.method = "fdr")##
## Pairwise comparisons using Wilcoxon rank sum test
##
## data: mOTUs_species5$maxbin.017 and mOTUs_species5$Zones
##
## Lam T10
## T10 0.333 -
## T4 0.095 0.548
##
## P value adjustment method: fdr
pairwise.wilcox.test(mOTUs_species5$`maxbin.019`, mOTUs_species5$Zones,p.adjust.method = "fdr")##
## Pairwise comparisons using Wilcoxon rank sum test
##
## data: mOTUs_species5$maxbin.019 and mOTUs_species5$Zones
##
## Lam T10
## T10 0.083 -
## T4 0.690 0.048
##
## P value adjustment method: fdr
pairwise.wilcox.test(mOTUs_species5$`maxbin.020`, mOTUs_species5$Zones,p.adjust.method = "fdr")##
## Pairwise comparisons using Wilcoxon rank sum test
##
## data: mOTUs_species5$maxbin.020 and mOTUs_species5$Zones
##
## Lam T10
## T10 0.048 -
## T4 0.024 0.056
##
## P value adjustment method: fdr
pairwise.wilcox.test(mOTUs_species5$`maxbin.026`, mOTUs_species5$Zones,p.adjust.method = "fdr")##
## Pairwise comparisons using Wilcoxon rank sum test
##
## data: mOTUs_species5$maxbin.026 and mOTUs_species5$Zones
##
## Lam T10
## T10 0.310 -
## T4 0.024 0.143
##
## P value adjustment method: fdr
pairwise.wilcox.test(mOTUs_species5$`maxbin.029`, mOTUs_species5$Zones,p.adjust.method = "fdr")##
## Pairwise comparisons using Wilcoxon rank sum test
##
## data: mOTUs_species5$maxbin.029 and mOTUs_species5$Zones
##
## Lam T10
## T10 0.548 -
## T4 0.048 0.083
##
## P value adjustment method: fdr
pairwise.wilcox.test(mOTUs_species5$`maxbin.033`, mOTUs_species5$Zones,p.adjust.method = "fdr")##
## Pairwise comparisons using Wilcoxon rank sum test
##
## data: mOTUs_species5$maxbin.033 and mOTUs_species5$Zones
##
## Lam T10
## T10 0.048 -
## T4 0.083 0.690
##
## P value adjustment method: fdr
pairwise.wilcox.test(mOTUs_species5$`maxbin.034`, mOTUs_species5$Zones,p.adjust.method = "fdr")##
## Pairwise comparisons using Wilcoxon rank sum test
##
## data: mOTUs_species5$maxbin.034 and mOTUs_species5$Zones
##
## Lam T10
## T10 0.226 -
## T4 0.024 0.310
##
## P value adjustment method: fdr
pairwise.wilcox.test(mOTUs_species5$`maxbin.043`, mOTUs_species5$Zones,p.adjust.method = "fdr")##
## Pairwise comparisons using Wilcoxon rank sum test
##
## data: mOTUs_species5$maxbin.043 and mOTUs_species5$Zones
##
## Lam T10
## T10 0.012 -
## T4 0.151 0.012
##
## P value adjustment method: fdr
pairwise.wilcox.test(mOTUs_species5$`maxbin.044`, mOTUs_species5$Zones,p.adjust.method = "fdr")##
## Pairwise comparisons using Wilcoxon rank sum test
##
## data: mOTUs_species5$maxbin.044 and mOTUs_species5$Zones
##
## Lam T10
## T10 0.048 -
## T4 0.083 0.151
##
## P value adjustment method: fdr
pairwise.wilcox.test(mOTUs_species5$`maxbin.056`, mOTUs_species5$Zones,p.adjust.method = "fdr")##
## Pairwise comparisons using Wilcoxon rank sum test
##
## data: mOTUs_species5$maxbin.056 and mOTUs_species5$Zones
##
## Lam T10
## T10 0.095 -
## T4 0.841 0.841
##
## P value adjustment method: fdr
pairwise.wilcox.test(mOTUs_species5$`maxbin.057`, mOTUs_species5$Zones,p.adjust.method = "fdr")##
## Pairwise comparisons using Wilcoxon rank sum test
##
## data: mOTUs_species5$maxbin.057 and mOTUs_species5$Zones
##
## Lam T10
## T10 1.000 -
## T4 0.048 0.226
##
## P value adjustment method: fdr
pairwise.wilcox.test(mOTUs_species5$`maxbin.062`, mOTUs_species5$Zones,p.adjust.method = "fdr")##
## Pairwise comparisons using Wilcoxon rank sum test
##
## data: mOTUs_species5$maxbin.062 and mOTUs_species5$Zones
##
## Lam T10
## T10 0.421 -
## T4 0.048 0.421
##
## P value adjustment method: fdr
pairwise.wilcox.test(mOTUs_species5$`maxbin.067`, mOTUs_species5$Zones,p.adjust.method = "fdr")## Warning in wilcox.test.default(xi, xj, paired = paired, ...): cannot
## compute exact p-value with ties
## Warning in wilcox.test.default(xi, xj, paired = paired, ...): cannot
## compute exact p-value with ties
## Warning in wilcox.test.default(xi, xj, paired = paired, ...): cannot
## compute exact p-value with ties
##
## Pairwise comparisons using Wilcoxon rank sum test
##
## data: mOTUs_species5$maxbin.067 and mOTUs_species5$Zones
##
## Lam T10
## T10 0.121 -
## T4 0.067 0.067
##
## P value adjustment method: fdr
pairwise.wilcox.test(mOTUs_species5$`maxbin.071`, mOTUs_species5$Zones,p.adjust.method = "fdr")##
## Pairwise comparisons using Wilcoxon rank sum test
##
## data: mOTUs_species5$maxbin.071 and mOTUs_species5$Zones
##
## Lam T10
## T10 0.143 -
## T4 0.048 0.421
##
## P value adjustment method: fdr
pairwise.wilcox.test(mOTUs_species5$`maxbin.087`, mOTUs_species5$Zones,p.adjust.method = "fdr")## Warning in wilcox.test.default(xi, xj, paired = paired, ...): cannot
## compute exact p-value with ties
##
## Pairwise comparisons using Wilcoxon rank sum test
##
## data: mOTUs_species5$maxbin.087 and mOTUs_species5$Zones
##
## Lam T10
## T10 0.012 -
## T4 0.012 0.116
##
## P value adjustment method: fdr
pairwise.wilcox.test(mOTUs_species5$`maxbin.123`, mOTUs_species5$Zones,p.adjust.method = "fdr")## Warning in wilcox.test.default(xi, xj, paired = paired, ...): cannot
## compute exact p-value with ties
##
## Pairwise comparisons using Wilcoxon rank sum test
##
## data: mOTUs_species5$maxbin.123 and mOTUs_species5$Zones
##
## Lam T10
## T10 0.016 -
## T4 0.016 0.016
##
## P value adjustment method: fdr
pairwise.wilcox.test(mOTUs_species5$`maxbin.124`, mOTUs_species5$Zones,p.adjust.method = "fdr")##
## Pairwise comparisons using Wilcoxon rank sum test
##
## data: mOTUs_species5$maxbin.124 and mOTUs_species5$Zones
##
## Lam T10
## T10 0.226 -
## T4 0.095 0.841
##
## P value adjustment method: fdr
pairwise.wilcox.test(mOTUs_species5$`maxbin.133`, mOTUs_species5$Zones,p.adjust.method = "fdr")##
## Pairwise comparisons using Wilcoxon rank sum test
##
## data: mOTUs_species5$maxbin.133 and mOTUs_species5$Zones
##
## Lam T10
## T10 0.690 -
## T4 0.095 0.143
##
## P value adjustment method: fdr
pairwise.wilcox.test(mOTUs_species5$`maxbin.135`, mOTUs_species5$Zones,p.adjust.method = "fdr")##
## Pairwise comparisons using Wilcoxon rank sum test
##
## data: mOTUs_species5$maxbin.135 and mOTUs_species5$Zones
##
## Lam T10
## T10 0.226 -
## T4 0.548 0.095
##
## P value adjustment method: fdr
pairwise.wilcox.test(mOTUs_species5$`maxbin.141`, mOTUs_species5$Zones,p.adjust.method = "fdr")##
## Pairwise comparisons using Wilcoxon rank sum test
##
## data: mOTUs_species5$maxbin.141 and mOTUs_species5$Zones
##
## Lam T10
## T10 0.024 -
## T4 0.226 0.690
##
## P value adjustment method: fdr
pairwise.wilcox.test(mOTUs_species5$`maxbin.151`, mOTUs_species5$Zones,p.adjust.method = "fdr")##
## Pairwise comparisons using Wilcoxon rank sum test
##
## data: mOTUs_species5$maxbin.151 and mOTUs_species5$Zones
##
## Lam T10
## T10 0.024 -
## T4 0.548 0.083
##
## P value adjustment method: fdr
pairwise.wilcox.test(mOTUs_species5$`maxbin.152`, mOTUs_species5$Zones,p.adjust.method = "fdr")## Warning in wilcox.test.default(xi, xj, paired = paired, ...): cannot
## compute exact p-value with ties
##
## Pairwise comparisons using Wilcoxon rank sum test
##
## data: mOTUs_species5$maxbin.152 and mOTUs_species5$Zones
##
## Lam T10
## T10 0.012 -
## T4 0.012 0.753
##
## P value adjustment method: fdr
pairwise.wilcox.test(mOTUs_species5$`maxbin.157`, mOTUs_species5$Zones,p.adjust.method = "fdr")## Warning in wilcox.test.default(xi, xj, paired = paired, ...): cannot
## compute exact p-value with ties
## Warning in wilcox.test.default(xi, xj, paired = paired, ...): cannot
## compute exact p-value with ties
##
## Pairwise comparisons using Wilcoxon rank sum test
##
## data: mOTUs_species5$maxbin.157 and mOTUs_species5$Zones
##
## Lam T10
## T10 0.024 -
## T4 0.094 0.086
##
## P value adjustment method: fdr
pairwise.wilcox.test(mOTUs_species5$`maxbin.166`, mOTUs_species5$Zones,p.adjust.method = "fdr")##
## Pairwise comparisons using Wilcoxon rank sum test
##
## data: mOTUs_species5$maxbin.166 and mOTUs_species5$Zones
##
## Lam T10
## T10 0.024 -
## T4 0.310 0.083
##
## P value adjustment method: fdr
Lactobacillus_johnsonii=ggplot(mOTUs_species5, aes(x=Zones, y=`Lactobacillus johnsonii`)) + geom_boxplot(fill=c("#999999", "#E69F00", "#56B4E9"))+geom_dotplot(binaxis='y', stackdir='center', dotsize=0.6)+
labs(title = "Lactobacillus johnsonii",x=NULL,y ="RPKM")+theme_classic()
Lactobacillus_johnsonii_boxplot=Lactobacillus_johnsonii+theme(plot.title = element_text(size = 24,hjust=0.5))+theme(axis.text.x = element_text(color="black",size = 20))+theme(axis.text.y = element_text(color="black",size = 20),axis.title.y=element_text(size=20))
Lactobacillus_johnsonii_boxplot+scale_y_continuous(expand = c(0,0), limits = c(0,45))+geom_segment(y=40,yend=40,x=1,xend=2)+geom_segment(y=42,yend=42,x=1,xend=3)+geom_text(x=2,y=42,label="**",size=10)+geom_text(x=1.5,y=40,label="*",size=10)+theme(plot.margin=unit(c(1,1,1.5,1.2),"cm"))+theme(axis.title.y = element_text(margin = margin(t = 0, r = 20, b = 0, l = 0)))## `stat_bindot()` using `bins = 30`. Pick better value with `binwidth`.
maxbin009=ggplot(mOTUs_species5, aes(x=Zones, y=`maxbin.009`)) + geom_boxplot(fill=c("#999999", "#E69F00", "#56B4E9"))+geom_dotplot(binaxis='y', stackdir='center', dotsize=0.6)+labs(title = "CAG-1031 maxbin.009",x=NULL,y ="RPKM")+theme_classic()
maxbin009_boxplot=maxbin009+theme(plot.title = element_text(size = 24,hjust=0.5))+theme(axis.text.x = element_text(color="black",size = 20))+theme(axis.text.y = element_text(color="black",size = 20),axis.title.y=element_text(size=20))
maxbin009_boxplot+scale_y_continuous(expand = c(0,0), limits = c(0,82))+geom_segment(y=78,yend=78,x=1,xend=3)+geom_text(x=2,y=78,label="**",size=10)+theme(plot.margin=unit(c(1,1,1.5,1.2),"cm"))+theme(axis.title.y = element_text(margin = margin(t = 0, r = 20, b = 0, l = 0)))## `stat_bindot()` using `bins = 30`. Pick better value with `binwidth`.
maxbin121=ggplot(mOTUs_species5, aes(x=Zones, y=`maxbin.121`)) + geom_boxplot(fill=c("#999999", "#E69F00", "#56B4E9"))+geom_dotplot(binaxis='y', stackdir='center', dotsize=0.6)+labs(title = "CAG-1031 maxbin.121",x=NULL,y ="RPKM")+theme_classic()
maxbin121_boxplot=maxbin121+theme(plot.title = element_text(size = 24,hjust=0.5))+theme(axis.text.x = element_text(color="black",size = 20))+theme(axis.text.y = element_text(color="black",size = 20),axis.title.y=element_text(size=20))
maxbin121_boxplot+scale_y_continuous(expand = c(0,0), limits = c(0,64))+geom_segment(y=60,yend=60,x=1,xend=3)+geom_text(x=2,y=60,label="**",size=10)+theme(plot.margin=unit(c(1,1,1.5,1.2),"cm"))+theme(axis.title.y = element_text(margin = margin(t = 0, r = 20, b = 0, l = 0)))## `stat_bindot()` using `bins = 30`. Pick better value with `binwidth`.
Weissella_cibaria=ggplot(mOTUs_species5, aes(x=Zones, y=`Weissella cibaria`)) + geom_boxplot(fill=c("#999999", "#E69F00", "#56B4E9"))+geom_dotplot(binaxis='y', stackdir='center', dotsize=0.6)+labs(title = "Weissella cibaria",x=NULL,y ="RPKM")+theme_classic()
Weissella_cibaria_boxplot=Weissella_cibaria+theme(plot.title = element_text(size = 24,hjust=0.5))+theme(axis.text.x = element_text(color="black",size = 20))+theme(axis.text.y = element_text(color="black",size = 20),axis.title.y=element_text(size=20))
Weissella_cibaria_boxplot+scale_y_continuous(expand = c(0,0), limits = c(0,41))+geom_segment(y=37,yend=37,x=2,xend=3)+geom_segment(y=39,yend=39,x=1,xend=3)+geom_text(x=2,y=39,label="*",size=12)+geom_text(x=2.5,y=37,label="*",size=10)+theme(plot.margin=unit(c(1,1,1.5,1.2),"cm"))+theme(axis.title.y = element_text(margin = margin(t = 0, r = 20, b = 0, l = 0)))## `stat_bindot()` using `bins = 30`. Pick better value with `binwidth`.
Lactococcus_lactis_A=ggplot(mOTUs_species5, aes(x=Zones, y=`Lactococcus lactis_A`)) + geom_boxplot(fill=c("#999999", "#E69F00", "#56B4E9"))+geom_dotplot(binaxis='y', stackdir='center', dotsize=0.6)+labs(title = "Lactococcus lactis_A",x=NULL,y ="RPKM")+theme_classic()
Lactococcus_lactis_A_boxplot=Lactococcus_lactis_A+theme(plot.title = element_text(size = 24,hjust=0.5))+theme(axis.text.x = element_text(color="black",size = 20))+theme(axis.text.y = element_text(color="black",size = 20),axis.title.y=element_text(size=20))
Lactococcus_lactis_A_boxplot+scale_y_continuous(expand = c(0,0), limits = c(0,40))+geom_segment(y=36,yend=36,x=2,xend=3)+geom_segment(y=38,yend=38,x=1,xend=3)+geom_text(x=2,y=38,label="**",size=10)+geom_text(x=2.5,y=36,label="**",size=10)+theme(plot.margin=unit(c(1,1,1.5,1.2),"cm"))+theme(axis.title.y = element_text(margin = margin(t = 0, r = 20, b = 0, l = 0)))## `stat_bindot()` using `bins = 30`. Pick better value with `binwidth`.
Bacteroides_thetaiotaomicron=ggplot(mOTUs_species5, aes(x=Zones, y=`Bacteroides thetaiotaomicron`)) + geom_boxplot(fill=c("#999999", "#E69F00", "#56B4E9"))+geom_dotplot(binaxis='y', stackdir='center', dotsize=0.6)+
labs(title = "Bacteroides thetaiotaomicron",x=NULL,y ="RPKM")+theme_classic()
Bacteroides_thetaiotaomicron_boxplot=Bacteroides_thetaiotaomicron+theme(plot.title = element_text(size = 24,hjust=0.5))+theme(axis.text.x = element_text(color="black",size = 20))+theme(axis.text.y = element_text(color="black",size = 20),axis.title.y=element_text(size=20))
Bacteroides_thetaiotaomicron_boxplot+scale_y_continuous(expand = c(0,0), limits = c(0,22))+geom_segment(y=19,yend=19,x=1,xend=3)+geom_text(x=2,y=20,label="p=0.055",size=6)+theme(plot.margin=unit(c(1,1,1.5,1.2),"cm"))+theme(axis.title.y = element_text(margin = margin(t = 0, r = 20, b = 0, l = 0)))## `stat_bindot()` using `bins = 30`. Pick better value with `binwidth`.
unique_PCs=read.csv("PCs_RPKM.csv",header=T,sep=",")
rownames(unique_PCs)=unique_PCs[,1]
unique_PCs=unique_PCs[,-1]
t_unique_PCs=t(unique_PCs)
# Applying a cube root transformation
d.clr_PCs = log2(t_unique_PCs+1)
# Calculating Bray-Curtis dissimilarity
Bray_PCs=vegdist(d.clr_PCs,"bray")
# Applying a non-constrained ordination (PCoA) on the dissimilarity matrix
unique_PCs_bray=capscale(Bray_PCs~-1)
# Extracting the percentage explained by the first two dimensions and automatically adding them to the axes titles
unique_PCs_bray_eig = eigenvals(unique_PCs_bray)
percentage_variance_explained <- unique_PCs_bray_eig/sum(unique_PCs_bray_eig)
sum_percentage_variance_explained <- cumsum(unique_PCs_bray_eig / sum(unique_PCs_bray_eig))
xlabel= as.numeric(format(round((percentage_variance_explained[1]*100), 2), nsmall = 2))
xlabel= sprintf("%.2f %%", xlabel)
xlabel= paste ("PCo1 (", xlabel, "of variation explained", ")")
ylabel= as.numeric(format(round((percentage_variance_explained[2]*100), 2), nsmall = 2))
ylabel= sprintf("%.2f %%", ylabel)
ylabel= paste ("PCo2 (", ylabel,"of variation explained", ")")
plot(unique_PCs_bray, type="n", xlab="", ylab="",ylim=c(-1, 1),xlim=c(-1,1),cex.axis=0.6, tck = -0.01, mgp = c(3, 1, 0), xaxp = c(-2, 2, 4), panel.first=grid(col = "white",lty=0))
title(ylab=ylabel, line=2, cex.lab=1)
title(xlab=xlabel, line=2, cex.lab=1)
abline(h=0, v=0, col = "white", lty = 1, lwd = 1.5)
abline(h=-10:10, v=-10:10, col = "lightgray", lty = "dotted", lwd = 0.5)
par(lty=2)
colvec_all <- c("#999999", "#999999","#999999","#999999","#999999","#E69F00","#E69F00","#E69F00","#E69F00","#E69F00", "#56B4E9","#56B4E9","#56B4E9","#56B4E9","#56B4E9")
points(unique_PCs_bray,cex= 2.5,pch=21, col="black", bg= colvec_all, lwd = 1)
# Adding a legend
legend(-1.5,1, pt.cex=2.5 , pt.lwd = 1,c("Lam", "T10", "T4"),bty = "n" ,pch = 21,col="black",pt.bg = c("#999999", "#E69F00", "#56B4E9"), cex = 1.5)
ordihull(unique_PCs_bray, Colvec_expedition$Zones, display = "sites", label = T)# ANOSIM statistics
PCs_anosim=anosim(Bray_PCs, permutations = 999, grouping=Colvec_expedition$Zones)
summary(PCs_anosim)##
## Call:
## anosim(x = Bray_PCs, grouping = Colvec_expedition$Zones, permutations = 999)
## Dissimilarity: bray
##
## ANOSIM statistic R: 0.4151
## Significance: 0.001
##
## Permutation: free
## Number of permutations: 999
##
## Upper quantiles of permutations (null model):
## 90% 95% 97.5% 99%
## 0.127 0.176 0.200 0.236
##
## Dissimilarity ranks between and within classes:
## 0% 25% 50% 75% 100% N
## Between 8 38.50 59.0 82.50 104 75
## Lam 1 3.50 9.0 16.25 26 10
## T10 4 9.25 29.5 74.25 95 10
## T4 16 42.00 65.0 72.50 105 10
#ANOSIM statistic R: 0.4151
#Significance: 0.001
# Each data point indicates an individual mouse sample.selected_functions=read.csv("selected_functions_MAGs.csv",header=T,row.names=NULL,sep=",")
row.names(selected_functions)=selected_functions[,1]
selected_functions=selected_functions[,-1]
#Lam and T4
selected_functions1 <- selected_functions[,c(idx_Lam, idx_T4)]
clin.sub_selected_functions1 <- data.frame(ID = 1:ncol(selected_functions1), group = rep(c("Lam","T4"), each=5))
fx <- function(x) {
res <-wilcox.test(x ~clin.sub_selected_functions1$group)
return(res$p.value)
}
selected_functions.wilcox1 <-apply(as.matrix(selected_functions1), 1, fx)
selected_functions.wilcox1 <- data.frame(selected_functions.wilcox1 )
colnames(selected_functions.wilcox1) <- c("p.value")
write.csv(selected_functions.wilcox1, file = "selected_functions.wilcox0.05_Lam_T4.csv")
rownames(selected_functions.wilcox1)[selected_functions.wilcox1$p.value< 0.01]## [1] "prephenate dehydratase"
## [2] "anthranilate synthase component II "
## [3] "indole-3-glycerol phosphate synthase "
## [4] "pyridoxine 4-dehydrogenase"
## [5] "dihydrofolate reductase"
## [6] "dihydroneopterin aldolase "
## [7] "lactocepin"
## [8] "lactose-specific components"
## [9] "N-acetylgalactosamine-specific components"
## [10] "galactosamine-specific components"
## [11] " glucitol/sorbitol-specific components"
## [12] "cellobiose phosphorylase"
rownames(selected_functions.wilcox1)[selected_functions.wilcox1$p.value< 0.05]## [1] "choloylglycine hydrolase"
## [2] "prephenate dehydratase"
## [3] "anthranilate synthase component II "
## [4] "indole-3-glycerol phosphate synthase "
## [5] "chorismate synthase"
## [6] "pyridoxine 4-dehydrogenase"
## [7] "dihydrofolate reductase"
## [8] "dihydroneopterin aldolase "
## [9] "lactocepin"
## [10] "lactose-specific components"
## [11] "N-acetylgalactosamine-specific components"
## [12] "galactosamine-specific components"
## [13] "galactitol-specific components"
## [14] "mannose-specific components "
## [15] "mannitol-specific components "
## [16] " glucitol/sorbitol-specific components"
## [17] "endo-1,4-beta-xylanase"
## [18] "arabinan endo-1,5-alpha-L-arabinosidase "
## [19] "mannan endo-1,4-beta-mannosidase"
## [20] "maltose phosphorylase"
## [21] "cellobiose phosphorylase"
#Lam and T10
selected_functions2 <- selected_functions[,c(idx_Lam, idx_T10)]
clin.sub_selected_functions2 <- data.frame(ID = 1:ncol(selected_functions2), group = rep(c("Lam","T10"), each=5))
fx <- function(x) {
res <-wilcox.test(x ~clin.sub_selected_functions2$group)
return(res$p.value)
}
selected_functions.wilcox2 <-apply(as.matrix(selected_functions2), 1, fx)
selected_functions.wilcox2 <- data.frame(selected_functions.wilcox2 )
colnames(selected_functions.wilcox2) <- c("p.value")
write.csv(selected_functions.wilcox2, file = "selected_functions.wilcox0.05_Lam_T10.csv")
rownames(selected_functions.wilcox2)[selected_functions.wilcox2$p.value< 0.01]## [1] "indole-3-glycerol phosphate synthase "
rownames(selected_functions.wilcox2)[selected_functions.wilcox2$p.value< 0.05]## [1] "prephenate dehydratase"
## [2] "indole-3-glycerol phosphate synthase "
## [3] "chorismate synthase"
## [4] "dihydroneopterin aldolase "
## [5] "lactocepin"
## [6] "lactose-specific components"
## [7] "N-acetylgalactosamine-specific components"
## [8] "galactosamine-specific components"
## [9] "mannose-specific components "
## [10] " glucitol/sorbitol-specific components"
## [11] "maltose phosphorylase"
## [12] "cellobiose phosphorylase"
#T10 and T4
selected_functions3 <-selected_functions[,c(idx_T4, idx_T10)]
clin.sub_selected_functions3 <- data.frame(ID = 1:ncol(selected_functions3), group = rep(c("T4","T10"), each=5))
fx <- function(x) {
res <-wilcox.test(x ~clin.sub_selected_functions3$group)
return(res$p.value)
}
selected_functions.wilcox3 <-apply(as.matrix(selected_functions3), 1, fx)
selected_functions.wilcox3 <- data.frame(selected_functions.wilcox3 )
colnames(selected_functions.wilcox3) <- c("p.value")
write.csv(selected_functions.wilcox3, file = "selected_functions.wilcox0.05_T4_T10.csv")
rownames(selected_functions.wilcox3)[selected_functions.wilcox3$p.value< 0.05]## character(0)
rownames(selected_functions.wilcox3)[selected_functions.wilcox3$p.value< 0.05]## character(0)
selected_functions=read.csv("selected_functions_MAGs.csv",header=T,row.names=NULL,sep=",")
row.names(selected_functions)=selected_functions[,1]
selected_functions=selected_functions[,-1]
selected_functions4=t(selected_functions)
selected_functions5=cbind(selected_functions4, Colvec_expedition)
pairwise.wilcox.test(as.numeric(selected_functions5$`choloylglycine hydrolase`),selected_functions5$Zones,p.adjust.method = "fdr")##
## Pairwise comparisons using Wilcoxon rank sum test
##
## data: as.numeric(selected_functions5$`choloylglycine hydrolase`) and selected_functions5$Zones
##
## Lam T10
## T10 0.083 -
## T4 0.083 0.690
##
## P value adjustment method: fdr
pairwise.wilcox.test(as.numeric(selected_functions5$`prephenate dehydratase`), selected_functions5$Zones,p.adjust.method = "fdr")##
## Pairwise comparisons using Wilcoxon rank sum test
##
## data: as.numeric(selected_functions5$`prephenate dehydratase`) and selected_functions5$Zones
##
## Lam T10
## T10 0.024 -
## T4 0.024 1.000
##
## P value adjustment method: fdr
pairwise.wilcox.test(as.numeric(selected_functions5$`anthranilate synthase component II `), mOTUs_species5$Zones,p.adjust.method = "fdr")##
## Pairwise comparisons using Wilcoxon rank sum test
##
## data: as.numeric(selected_functions5$`anthranilate synthase component II `) and mOTUs_species5$Zones
##
## Lam T10
## T10 0.690 -
## T4 0.024 0.631
##
## P value adjustment method: fdr
pairwise.wilcox.test(as.numeric(selected_functions5$`indole-3-glycerol phosphate synthase `), selected_functions5$Zones,p.adjust.method = "fdr")##
## Pairwise comparisons using Wilcoxon rank sum test
##
## data: as.numeric(selected_functions5$`indole-3-glycerol phosphate synthase `) and selected_functions5$Zones
##
## Lam T10
## T10 0.012 -
## T4 0.012 1.000
##
## P value adjustment method: fdr
pairwise.wilcox.test(as.numeric(selected_functions5$`chorismate synthase`), selected_functions5$Zones,p.adjust.method = "fdr")##
## Pairwise comparisons using Wilcoxon rank sum test
##
## data: as.numeric(selected_functions5$`chorismate synthase`) and selected_functions5$Zones
##
## Lam T10
## T10 0.048 -
## T4 0.048 0.222
##
## P value adjustment method: fdr
pairwise.wilcox.test(as.numeric(selected_functions5$`pyridoxine 4-dehydrogenase`),selected_functions5$Zones,p.adjust.method = "fdr")##
## Pairwise comparisons using Wilcoxon rank sum test
##
## data: as.numeric(selected_functions5$`pyridoxine 4-dehydrogenase`) and selected_functions5$Zones
##
## Lam T10
## T10 0.083 -
## T4 0.024 0.548
##
## P value adjustment method: fdr
pairwise.wilcox.test(as.numeric(selected_functions5$`dihydrofolate reductase`), selected_functions5$Zones,p.adjust.method = "fdr")##
## Pairwise comparisons using Wilcoxon rank sum test
##
## data: as.numeric(selected_functions5$`dihydrofolate reductase`) and selected_functions5$Zones
##
## Lam T10
## T10 0.143 -
## T4 0.024 0.151
##
## P value adjustment method: fdr
pairwise.wilcox.test(as.numeric(selected_functions5$`dihydroneopterin aldolase `),selected_functions5$Zones,p.adjust.method = "fdr")##
## Pairwise comparisons using Wilcoxon rank sum test
##
## data: as.numeric(selected_functions5$`dihydroneopterin aldolase `) and selected_functions5$Zones
##
## Lam T10
## T10 0.024 -
## T4 0.024 1.000
##
## P value adjustment method: fdr
pairwise.wilcox.test(as.numeric(selected_functions5$lactocepin), selected_functions5$Zones,p.adjust.method = "fdr")##
## Pairwise comparisons using Wilcoxon rank sum test
##
## data: as.numeric(selected_functions5$lactocepin) and selected_functions5$Zones
##
## Lam T10
## T10 0.024 -
## T4 0.024 0.548
##
## P value adjustment method: fdr
pairwise.wilcox.test(as.numeric(selected_functions5$`lactose-specific components`), selected_functions5$Zones,p.adjust.method = "fdr")##
## Pairwise comparisons using Wilcoxon rank sum test
##
## data: as.numeric(selected_functions5$`lactose-specific components`) and selected_functions5$Zones
##
## Lam T10
## T10 0.024 -
## T4 0.024 0.548
##
## P value adjustment method: fdr
pairwise.wilcox.test(as.numeric(selected_functions5$`N-acetylgalactosamine-specific components`), selected_functions5$Zones,p.adjust.method = "fdr")##
## Pairwise comparisons using Wilcoxon rank sum test
##
## data: as.numeric(selected_functions5$`N-acetylgalactosamine-specific components`) and selected_functions5$Zones
##
## Lam T10
## T10 0.024 -
## T4 0.024 0.548
##
## P value adjustment method: fdr
pairwise.wilcox.test(as.numeric(selected_functions5$`galactosamine-specific components`), selected_functions5$Zones,p.adjust.method = "fdr")##
## Pairwise comparisons using Wilcoxon rank sum test
##
## data: as.numeric(selected_functions5$`galactosamine-specific components`) and selected_functions5$Zones
##
## Lam T10
## T10 0.024 -
## T4 0.024 1.000
##
## P value adjustment method: fdr
pairwise.wilcox.test(as.numeric(selected_functions5$`galactitol-specific components`), selected_functions5$Zones,p.adjust.method = "fdr")##
## Pairwise comparisons using Wilcoxon rank sum test
##
## data: as.numeric(selected_functions5$`galactitol-specific components`) and selected_functions5$Zones
##
## Lam T10
## T10 0.143 -
## T4 0.048 0.548
##
## P value adjustment method: fdr
pairwise.wilcox.test(as.numeric(selected_functions5$`mannose-specific components `), selected_functions5$Zones,p.adjust.method = "BH")##
## Pairwise comparisons using Wilcoxon rank sum test
##
## data: as.numeric(selected_functions5$`mannose-specific components `) and selected_functions5$Zones
##
## Lam T10
## T10 0.048 -
## T4 0.048 0.151
##
## P value adjustment method: BH
pairwise.wilcox.test(as.numeric(selected_functions5$`mannitol-specific components `), selected_functions5$Zones,p.adjust.method = "fdr")##
## Pairwise comparisons using Wilcoxon rank sum test
##
## data: as.numeric(selected_functions5$`mannitol-specific components `) and selected_functions5$Zones
##
## Lam T10
## T10 0.222 -
## T4 0.095 0.222
##
## P value adjustment method: fdr
pairwise.wilcox.test(as.numeric(selected_functions5$` glucitol/sorbitol-specific components`), selected_functions5$Zones,p.adjust.method = "fdr")##
## Pairwise comparisons using Wilcoxon rank sum test
##
## data: as.numeric(selected_functions5$` glucitol/sorbitol-specific components`) and selected_functions5$Zones
##
## Lam T10
## T10 0.048 -
## T4 0.024 0.841
##
## P value adjustment method: fdr
pairwise.wilcox.test(as.numeric(selected_functions5$`endo-1,4-beta-xylanase`), selected_functions5$Zones,p.adjust.method = "fdr")##
## Pairwise comparisons using Wilcoxon rank sum test
##
## data: as.numeric(selected_functions5$`endo-1,4-beta-xylanase`) and selected_functions5$Zones
##
## Lam T10
## T10 0.143 -
## T4 0.048 0.841
##
## P value adjustment method: fdr
pairwise.wilcox.test(as.numeric(selected_functions5$`arabinan endo-1,5-alpha-L-arabinosidase `), selected_functions3$Zones,p.adjust.method = "fdr")##
## Pairwise comparisons using Wilcoxon rank sum test
##
## data: as.numeric(selected_functions5$`arabinan endo-1,5-alpha-L-arabinosidase `) and selected_functions3$Zones
##
## <0 x 0 matrix>
##
## P value adjustment method: fdr
pairwise.wilcox.test(as.numeric(selected_functions5$`mannan endo-1,4-beta-mannosidase`), selected_functions5$Zones,p.adjust.method = "fdr")##
## Pairwise comparisons using Wilcoxon rank sum test
##
## data: as.numeric(selected_functions5$`mannan endo-1,4-beta-mannosidase`) and selected_functions5$Zones
##
## Lam T10
## T10 0.143 -
## T4 0.095 0.841
##
## P value adjustment method: fdr
pairwise.wilcox.test(as.numeric(selected_functions5$`maltose phosphorylase`), selected_functions5$Zones,p.adjust.method = "fdr")##
## Pairwise comparisons using Wilcoxon rank sum test
##
## data: as.numeric(selected_functions5$`maltose phosphorylase`) and selected_functions5$Zones
##
## Lam T10
## T10 0.048 -
## T4 0.048 0.095
##
## P value adjustment method: fdr
pairwise.wilcox.test(as.numeric(selected_functions5$`cellobiose phosphorylase`), selected_functions5$Zones,p.adjust.method = "fdr")##
## Pairwise comparisons using Wilcoxon rank sum test
##
## data: as.numeric(selected_functions5$`cellobiose phosphorylase`) and selected_functions5$Zones
##
## Lam T10
## T10 0.048 -
## T4 0.024 0.421
##
## P value adjustment method: fdr
cal_z_score <- function(x){
(x - mean(x)) / sd(x)
}
selected_functions_norm <- t(apply(selected_functions, 1, cal_z_score))
my_sample_col <- data.frame(sample = rep(c("Lam","T10","T4"), each=5))
row.names(my_sample_col) <- colnames(selected_functions)
pheatmap(selected_functions_norm,cluster_rows=FALSE,cluster_cols=FALSE,border_color = NA,annotation_col = my_sample_col)#B.Differential abundances of bacterial genera (p<0.05 by Wilcoxon signed-rank test) in either two groups are indicated in red. Each row representing a unique bacterial genus was Z-score normalized. Bacterial genera on the y-axis are clustered using Euclidean distances.
#use pheatmap
cal_z_score <- function(x){
(x - mean(x)) / sd(x)
}
genus_singleM_norm <- t(apply(genus_singleM, 1, cal_z_score))
pheatmap(genus_singleM_norm,cluster_cols=FALSE, legend = TRUE,annotation_col = my_sample_col)singleM_genus_rare=data_frame(singleM_genus$Clostridium_M,singleM_genus$Eubacterium_I,singleM_genus$Slackia,singleM_genus$`CAG-791`,singleM_genus$Pseudooceanicola,singleM_genus$Ruminococcus,singleM_genus$Flavonifractor,singleM_genus$Zones)
rownames(singleM_genus_rare)=rownames(singleM_genus)## Warning: Setting row names on a tibble is deprecated.
colnames(singleM_genus_rare)=c("Clostridium_M","Eubacterium_I","Slackia","CAG-791","Pseudooceanicola","Ruminococcus","Flavonifractor","Zones")
singleM_genus_rare1=gather(singleM_genus_rare,key="Rare",value = "abundances",-Zones)
singleM_genus_rare1$Rare=factor(singleM_genus_rare1$Rare,levels=c("Clostridium_M","Eubacterium_I","Flavonifractor","Slackia","CAG-791","Ruminococcus","Pseudooceanicola"))
genus_rare_boxplot=ggplot(singleM_genus_rare1, aes(x=Rare, y=abundances,fill=factor(Zones)))+theme_classic()+scale_fill_manual(values=c("#999999", "#E69F00", "#56B4E9"))+geom_boxplot(position=position_dodge(0.8))+ geom_dotplot(binaxis='y', stackdir='center', position=position_dodge(0.8),binpositions = "bygroup", dotsize=0.4)+labs(title =NULL,x=NULL,y = "Class Clostridia \n relative abundances(%)")+theme_classic()
genus_rare_boxplot=genus_rare_boxplot+theme(plot.title = element_text(size = 28,hjust=0.5))+theme(axis.text.x = element_text(color="black",size = 10,angle=45,hjust = 1))+theme(axis.text.y = element_text(color="black",size = 10),axis.title.y=element_text(size=12))
genus_rare_boxplot+scale_y_continuous(expand = c(0,0), limits = c(0,0.6))+geom_segment(y=0.5,yend=0.5,x=0.75,xend=1)+geom_text(x=0.875,y=0.5,label="*",size=6)+geom_segment(y=0.29,yend=0.29,x=1,xend=1.25)+geom_text(x=1.125,y=0.29,label="*",size=6)+geom_segment(y=0.31,yend=0.31,x=2,xend=2.25)+geom_text(x=2.125,y=0.31,label="**",size=6)+geom_segment(y=0.23,yend=0.23,x=2.75,xend=3)+geom_text(x=2.875,y=0.23,label="*",size=6)+geom_segment(y=0.11,yend=0.11,x=3.75,xend=4)+geom_text(x=3.875,y=0.11,label="*",size=6)+geom_segment(y=0.11,yend=0.11,x=5,xend=5.25)+geom_text(x=5.125,y=0.11,label="**",size=6)+geom_segment(y=0.135,yend=0.135,x=5.75,xend=6)+geom_text(x=5.875,y=0.135,label="*",size=6)+geom_segment(y=0.03,yend=0.03,x=6.75,xend=7)+geom_text(x=6.875,y=0.03,label="*",size=6)+theme(axis.title.y = element_text(margin = margin(t = 0, r = 20, b = 0, l = 0)))## `stat_bindot()` using `bins = 30`. Pick better value with `binwidth`.
#Differential abundances of bacterial species (p<0.05 by Wilcoxon signed-rank test) in either two groups are indicated in red. Each row representing a unique bacterial genus was Z-score normalized.
mOTUs_species=read.table("MAGs_species.txt",header=T,row.names=1,sep="\t")
mOTUs_species$Maxbins=as.factor(rownames((mOTUs_species)))
mOTUs_species1=merge(mOTUs_species,mOTUs,by="Maxbins")
mOTUs_species2=mOTUs_species1[,-1]
row.names(mOTUs_species2)=mOTUs_species2[,1]
mOTUs_species3=mOTUs_species2[,-1]
mOTUs_species_norm <- t(apply(mOTUs_species3, 1, cal_z_score))
my_sample_col <- data.frame(sample = rep(c("Lam","T10","T4"), each=5))
row.names(my_sample_col) <- colnames(mOTUs_species3)
pheatmap(mOTUs_species_norm,cluster_cols=FALSE, annotation_col = my_sample_col,legend = TRUE)KEGG_carbohydrate=read.csv("mean_KEGGFUN_wilcox_0.05_carbohydrate_metabolism.csv",header=T,sep=",")
rownames(KEGG_carbohydrate)=KEGG_carbohydrate[,1]
KEGG_carbohydrate1=KEGG_carbohydrate[,-(1:2)]
cal_z_score <- function(x){
(x - mean(x)) / sd(x)
}
KEGG_carbohydrate1_norm <- t(apply(KEGG_carbohydrate1, 1, cal_z_score))
my_sample_col <- data.frame(sample = rep(c("Lam","T10","T4"), each=5))
row.names(my_sample_col) <- colnames(KEGG_carbohydrate1_norm)
pheatmap(KEGG_carbohydrate1_norm,annotation_col=my_sample_col,show_colnames = F, show_rownames = T,legend = TRUE, cluster_rows = T,cluster_cols = F,border_color = NA,cellwidth = 12, cellheight = 15, main="Carbohydrate Metabolism")#Lam and T4
KEGG_carbohydrate_functions1 <- KEGG_carbohydrate1[,c(idx_Lam, idx_T4)]
clin.sub_KEGG_carbohydrate_functions1 <- data.frame(ID = 1:ncol(KEGG_carbohydrate_functions1), group = rep(c("Lam","T4"), each=5))
fx <- function(x) {
res <-wilcox.test(x ~clin.sub_KEGG_carbohydrate_functions1$group)
return(res$p.value)
}
KEGG_carbohydrate_functions.wilcox1 <-apply(as.matrix(KEGG_carbohydrate_functions1), 1, fx)## Warning in wilcox.test.default(x = structure(c(0.617, 0, 0.228, 0,
## 0), .Names = c("Lam_4.x", : cannot compute exact p-value with ties
## Warning in wilcox.test.default(x = structure(c(0, 0, 0, 0, 0), .Names =
## c("Lam_4.x", : cannot compute exact p-value with ties
KEGG_carbohydrate_functions.wilcox1 <- data.frame(KEGG_carbohydrate_functions.wilcox1 )
colnames(KEGG_carbohydrate_functions.wilcox1) <- c("p.value")
write.csv(KEGG_carbohydrate_functions.wilcox1, file = "KEGG_carbohydrate_functions.wilcox0.05_Lam_T4.csv")
#Lam and T10
KEGG_carbohydrate_functions2 <- KEGG_carbohydrate1[,c(idx_Lam, idx_T10)]
clin.sub_KEGG_carbohydrate_functions2 <- data.frame(ID = 1:ncol(KEGG_carbohydrate_functions2), group = rep(c("Lam","T10"), each=5))
fx <- function(x) {
res <-wilcox.test(x ~clin.sub_KEGG_carbohydrate_functions2$group)
return(res$p.value)
}
KEGG_carbohydrate_functions.wilcox2 <-apply(as.matrix(KEGG_carbohydrate_functions2), 1, fx)## Warning in wilcox.test.default(x = structure(c(0.617, 0, 0.228, 0,
## 0), .Names = c("Lam_4.x", : cannot compute exact p-value with ties
## Warning in wilcox.test.default(x = structure(c(0, 0, 0, 0, 0), .Names =
## c("Lam_4.x", : cannot compute exact p-value with ties
KEGG_carbohydrate_functions.wilcox2 <- data.frame(KEGG_carbohydrate_functions.wilcox2 )
colnames(KEGG_carbohydrate_functions.wilcox2) <- c("p.value")
write.csv(KEGG_carbohydrate_functions.wilcox2, file = "KEGG_carbohydrate_functions.wilcox0.05_Lam_T10.csv")
#T10 and T4
KEGG_carbohydrate_functions3 <-KEGG_carbohydrate1[,c(idx_T4, idx_T10)]
clin.sub_KEGG_carbohydrate_functions3 <- data.frame(ID = 1:ncol(KEGG_carbohydrate_functions3), group = rep(c("T4","T10"), each=5))
fx <- function(x) {
res <-wilcox.test(x ~clin.sub_KEGG_carbohydrate_functions3$group)
return(res$p.value)
}
KEGG_carbohydrate_functions.wilcox3 <-apply(as.matrix(KEGG_carbohydrate_functions3), 1, fx)## Warning in wilcox.test.default(x = structure(c(0, 0.101, 1.04, 0.502,
## 0), .Names = c("T10_12.x", : cannot compute exact p-value with ties
## Warning in wilcox.test.default(x = structure(c(0, 0, 0, 0, 0.299), .Names =
## c("T10_12.x", : cannot compute exact p-value with ties
KEGG_carbohydrate_functions.wilcox3 <- data.frame(KEGG_carbohydrate_functions.wilcox3 )
colnames(KEGG_carbohydrate_functions.wilcox3) <- c("p.value")
write.csv(KEGG_carbohydrate_functions.wilcox3, file = "KEGG_carbohydrate_functions.wilcox0.05_T4_T10.csv")#Wilcoxon Rank Sum test with false discovery rates
KEGG_carbohydrate2=t(KEGG_carbohydrate1)
KEGG_carbohydrate3=cbind(KEGG_carbohydrate2, Colvec_expedition)
pairwise.wilcox.test(KEGG_carbohydrate3$`chitinase [EC:3.2.1.14]`, KEGG_carbohydrate3$Zones,p.adjust.method = "fdr")##
## Pairwise comparisons using Wilcoxon rank sum test
##
## data: KEGG_carbohydrate3$`chitinase [EC:3.2.1.14]` and KEGG_carbohydrate3$Zones
##
## Lam T10
## T10 0.143 -
## T4 0.048 0.421
##
## P value adjustment method: fdr
pairwise.wilcox.test(KEGG_carbohydrate3$`N-acetylmuramic acid 6-phosphate etherase [EC:4.2.-.-]`, KEGG_carbohydrate3$Zones,p.adjust.method = "fdr")##
## Pairwise comparisons using Wilcoxon rank sum test
##
## data: KEGG_carbohydrate3$`N-acetylmuramic acid 6-phosphate etherase [EC:4.2.-.-]` and KEGG_carbohydrate3$Zones
##
## Lam T10
## T10 0.048 -
## T4 0.048 1.000
##
## P value adjustment method: fdr
pairwise.wilcox.test(KEGG_carbohydrate3$`N-acetylneuraminate synthase [EC:2.5.1.56]`, KEGG_carbohydrate3$Zones,p.adjust.method = "fdr")##
## Pairwise comparisons using Wilcoxon rank sum test
##
## data: KEGG_carbohydrate3$`N-acetylneuraminate synthase [EC:2.5.1.56]` and KEGG_carbohydrate3$Zones
##
## Lam T10
## T10 0.690 -
## T4 0.690 0.095
##
## P value adjustment method: fdr
pairwise.wilcox.test(KEGG_carbohydrate3$`N-acylglucosamine 2-epimerase [EC:5.1.3.8]`, KEGG_carbohydrate3$Zones,p.adjust.method = "fdr")##
## Pairwise comparisons using Wilcoxon rank sum test
##
## data: KEGG_carbohydrate3$`N-acylglucosamine 2-epimerase [EC:5.1.3.8]` and KEGG_carbohydrate3$Zones
##
## Lam T10
## T10 0.226 -
## T4 0.024 1.000
##
## P value adjustment method: fdr
pairwise.wilcox.test(KEGG_carbohydrate3$`NAD-dependent deacetylase [EC:3.5.1.-]`, KEGG_carbohydrate3$Zones,p.adjust.method = "fdr")##
## Pairwise comparisons using Wilcoxon rank sum test
##
## data: KEGG_carbohydrate3$`NAD-dependent deacetylase [EC:3.5.1.-]` and KEGG_carbohydrate3$Zones
##
## Lam T10
## T10 0.056 -
## T4 0.024 0.056
##
## P value adjustment method: fdr
pairwise.wilcox.test(KEGG_carbohydrate3$`UDP-N-acetylmuramate dehydrogenase [EC:1.1.1.158]`, KEGG_carbohydrate3$Zones,p.adjust.method = "fdr")##
## Pairwise comparisons using Wilcoxon rank sum test
##
## data: KEGG_carbohydrate3$`UDP-N-acetylmuramate dehydrogenase [EC:1.1.1.158]` and KEGG_carbohydrate3$Zones
##
## Lam T10
## T10 0.048 -
## T4 0.048 0.690
##
## P value adjustment method: fdr
pairwise.wilcox.test(KEGG_carbohydrate3$`acetoin dehydrogenase [EC:1.1.1.5]`, KEGG_carbohydrate3$Zones,p.adjust.method = "fdr")## Warning in wilcox.test.default(xi, xj, paired = paired, ...): cannot
## compute exact p-value with ties
## Warning in wilcox.test.default(xi, xj, paired = paired, ...): cannot
## compute exact p-value with ties
## Warning in wilcox.test.default(xi, xj, paired = paired, ...): cannot
## compute exact p-value with ties
##
## Pairwise comparisons using Wilcoxon rank sum test
##
## data: KEGG_carbohydrate3$`acetoin dehydrogenase [EC:1.1.1.5]` and KEGG_carbohydrate3$Zones
##
## Lam T10
## T10 0.656 -
## T4 0.089 0.089
##
## P value adjustment method: fdr
pairwise.wilcox.test(KEGG_carbohydrate3$`succinate-semialdehyde dehydrogenase (NADP+) [EC:1.2.1.16]`, KEGG_carbohydrate3$Zones,p.adjust.method = "fdr")##
## Pairwise comparisons using Wilcoxon rank sum test
##
## data: KEGG_carbohydrate3$`succinate-semialdehyde dehydrogenase (NADP+) [EC:1.2.1.16]` and KEGG_carbohydrate3$Zones
##
## Lam T10
## T10 0.048 -
## T4 0.024 0.690
##
## P value adjustment method: fdr
pairwise.wilcox.test(KEGG_carbohydrate3$`hydroxymethylglutaryl-CoA synthase [EC:2.3.3.10]`, KEGG_carbohydrate3$Zones,p.adjust.method = "fdr")##
## Pairwise comparisons using Wilcoxon rank sum test
##
## data: KEGG_carbohydrate3$`hydroxymethylglutaryl-CoA synthase [EC:2.3.3.10]` and KEGG_carbohydrate3$Zones
##
## Lam T10
## T10 0.048 -
## T4 0.048 0.151
##
## P value adjustment method: fdr
pairwise.wilcox.test(KEGG_carbohydrate3$`isocitrate dehydrogenase (NAD+) [EC:1.1.1.41]`, KEGG_carbohydrate3$Zones,p.adjust.method = "fdr")##
## Pairwise comparisons using Wilcoxon rank sum test
##
## data: KEGG_carbohydrate3$`isocitrate dehydrogenase (NAD+) [EC:1.1.1.41]` and KEGG_carbohydrate3$Zones
##
## Lam T10
## T10 0.083 -
## T4 0.083 1.000
##
## P value adjustment method: fdr
pairwise.wilcox.test(KEGG_carbohydrate3$`citrate lyase subunit beta [EC:4.1.3.6]`, KEGG_carbohydrate3$Zones,p.adjust.method = "fdr")##
## Pairwise comparisons using Wilcoxon rank sum test
##
## data: KEGG_carbohydrate3$`citrate lyase subunit beta [EC:4.1.3.6]` and KEGG_carbohydrate3$Zones
##
## Lam T10
## T10 0.310 -
## T4 0.095 0.143
##
## P value adjustment method: fdr
pairwise.wilcox.test(KEGG_carbohydrate3$`succinyl-CoA synthetase alpha subunit [EC:6.2.1.5]`, KEGG_carbohydrate3$Zones,p.adjust.method = "fdr")##
## Pairwise comparisons using Wilcoxon rank sum test
##
## data: KEGG_carbohydrate3$`succinyl-CoA synthetase alpha subunit [EC:6.2.1.5]` and KEGG_carbohydrate3$Zones
##
## Lam T10
## T10 0.333 -
## T4 0.048 0.421
##
## P value adjustment method: fdr
pairwise.wilcox.test(KEGG_carbohydrate3$`L-iditol 2-dehydrogenase [EC:1.1.1.14]`, KEGG_carbohydrate3$Zones,p.adjust.method = "fdr")##
## Pairwise comparisons using Wilcoxon rank sum test
##
## data: KEGG_carbohydrate3$`L-iditol 2-dehydrogenase [EC:1.1.1.14]` and KEGG_carbohydrate3$Zones
##
## Lam T10
## T10 0.690 -
## T4 0.095 0.143
##
## P value adjustment method: fdr
pairwise.wilcox.test(KEGG_carbohydrate3$`L-rhamnose isomerase [EC:5.3.1.14]`, KEGG_carbohydrate3$Zones,p.adjust.method = "fdr")##
## Pairwise comparisons using Wilcoxon rank sum test
##
## data: KEGG_carbohydrate3$`L-rhamnose isomerase [EC:5.3.1.14]` and KEGG_carbohydrate3$Zones
##
## Lam T10
## T10 0.226 -
## T4 0.095 0.690
##
## P value adjustment method: fdr
pairwise.wilcox.test(KEGG_carbohydrate3$`mannan endo-1,4-beta-mannosidase [EC:3.2.1.78]`, KEGG_carbohydrate3$Zones,p.adjust.method = "fdr")##
## Pairwise comparisons using Wilcoxon rank sum test
##
## data: KEGG_carbohydrate3$`mannan endo-1,4-beta-mannosidase [EC:3.2.1.78]` and KEGG_carbohydrate3$Zones
##
## Lam T10
## T10 0.143 -
## T4 0.095 0.841
##
## P value adjustment method: fdr
pairwise.wilcox.test(KEGG_carbohydrate3$`PTS system, glucitol/sorbitol-specific IIA component [EC:2.7.1.69]`, KEGG_carbohydrate3$Zones,p.adjust.method = "fdr")##
## Pairwise comparisons using Wilcoxon rank sum test
##
## data: KEGG_carbohydrate3$`PTS system, glucitol/sorbitol-specific IIA component [EC:2.7.1.69]` and KEGG_carbohydrate3$Zones
##
## Lam T10
## T10 0.048 -
## T4 0.024 0.841
##
## P value adjustment method: fdr
pairwise.wilcox.test(KEGG_carbohydrate3$`PTS system, mannitol-specific IIA component [EC:2.7.1.69]`, KEGG_carbohydrate3$Zones,p.adjust.method = "fdr")##
## Pairwise comparisons using Wilcoxon rank sum test
##
## data: KEGG_carbohydrate3$`PTS system, mannitol-specific IIA component [EC:2.7.1.69]` and KEGG_carbohydrate3$Zones
##
## Lam T10
## T10 0.222 -
## T4 0.095 0.222
##
## P value adjustment method: fdr
pairwise.wilcox.test(KEGG_carbohydrate3$`mannose-6-phosphate isomerase [EC:5.3.1.8]`, KEGG_carbohydrate3$Zones,p.adjust.method = "fdr")##
## Pairwise comparisons using Wilcoxon rank sum test
##
## data: KEGG_carbohydrate3$`mannose-6-phosphate isomerase [EC:5.3.1.8]` and KEGG_carbohydrate3$Zones
##
## Lam T10
## T10 0.226 -
## T4 0.095 0.548
##
## P value adjustment method: fdr
pairwise.wilcox.test(KEGG_carbohydrate3$`PTS system, mannose-specific IIA component [EC:2.7.1.69]`, KEGG_carbohydrate3$Zones,p.adjust.method = "fdr")##
## Pairwise comparisons using Wilcoxon rank sum test
##
## data: KEGG_carbohydrate3$`PTS system, mannose-specific IIA component [EC:2.7.1.69]` and KEGG_carbohydrate3$Zones
##
## Lam T10
## T10 0.024 -
## T4 0.024 0.548
##
## P value adjustment method: fdr
pairwise.wilcox.test(KEGG_carbohydrate3$`PTS system, mannose-specific IIB component [EC:2.7.1.69]`, KEGG_carbohydrate3$Zones,p.adjust.method = "fdr")##
## Pairwise comparisons using Wilcoxon rank sum test
##
## data: KEGG_carbohydrate3$`PTS system, mannose-specific IIB component [EC:2.7.1.69]` and KEGG_carbohydrate3$Zones
##
## Lam T10
## T10 0.048 -
## T4 0.048 0.222
##
## P value adjustment method: fdr
pairwise.wilcox.test(KEGG_carbohydrate3$`PTS system, mannose-specific IIC component`, KEGG_carbohydrate3$Zones,p.adjust.method = "fdr")##
## Pairwise comparisons using Wilcoxon rank sum test
##
## data: KEGG_carbohydrate3$`PTS system, mannose-specific IIC component` and KEGG_carbohydrate3$Zones
##
## Lam T10
## T10 0.048 -
## T4 0.048 0.151
##
## P value adjustment method: fdr
pairwise.wilcox.test(KEGG_carbohydrate3$`galactose-6-phosphate isomerase [EC:5.3.1.26]`, KEGG_carbohydrate3$Zones,p.adjust.method = "fdr")##
## Pairwise comparisons using Wilcoxon rank sum test
##
## data: KEGG_carbohydrate3$`galactose-6-phosphate isomerase [EC:5.3.1.26]` and KEGG_carbohydrate3$Zones
##
## Lam T10
## T10 0.048 -
## T4 0.048 0.690
##
## P value adjustment method: fdr
pairwise.wilcox.test(KEGG_carbohydrate3$`PTS system, galactitol-specific IIA component [EC:2.7.1.69]`,KEGG_carbohydrate3$Zones,p.adjust.method = "fdr")##
## Pairwise comparisons using Wilcoxon rank sum test
##
## data: KEGG_carbohydrate3$`PTS system, galactitol-specific IIA component [EC:2.7.1.69]` and KEGG_carbohydrate3$Zones
##
## Lam T10
## T10 0.024 -
## T4 0.024 0.421
##
## P value adjustment method: fdr
pairwise.wilcox.test(KEGG_carbohydrate3$`PTS system, galactitol-specific IIB component [EC:2.7.1.69]`, KEGG_carbohydrate3$Zones,p.adjust.method = "fdr")##
## Pairwise comparisons using Wilcoxon rank sum test
##
## data: KEGG_carbohydrate3$`PTS system, galactitol-specific IIB component [EC:2.7.1.69]` and KEGG_carbohydrate3$Zones
##
## Lam T10
## T10 0.048 -
## T4 0.048 0.310
##
## P value adjustment method: fdr
pairwise.wilcox.test(KEGG_carbohydrate3$`PTS system, galactitol-specific IIC component`, KEGG_carbohydrate3$Zones,p.adjust.method = "fdr")##
## Pairwise comparisons using Wilcoxon rank sum test
##
## data: KEGG_carbohydrate3$`PTS system, galactitol-specific IIC component` and KEGG_carbohydrate3$Zones
##
## Lam T10
## T10 0.222 -
## T4 0.083 0.083
##
## P value adjustment method: fdr
pairwise.wilcox.test(KEGG_carbohydrate3$`PTS system, galactosamine-specific IIB component [EC:2.7.1.69]`, KEGG_carbohydrate3$Zones,p.adjust.method = "fdr")##
## Pairwise comparisons using Wilcoxon rank sum test
##
## data: KEGG_carbohydrate3$`PTS system, galactosamine-specific IIB component [EC:2.7.1.69]` and KEGG_carbohydrate3$Zones
##
## Lam T10
## T10 0.024 -
## T4 0.024 0.841
##
## P value adjustment method: fdr
pairwise.wilcox.test(KEGG_carbohydrate3$`PTS system, galactosamine-specific IIC component`, KEGG_carbohydrate3$Zones,p.adjust.method = "fdr")##
## Pairwise comparisons using Wilcoxon rank sum test
##
## data: KEGG_carbohydrate3$`PTS system, galactosamine-specific IIC component` and KEGG_carbohydrate3$Zones
##
## Lam T10
## T10 0.024 -
## T4 0.024 0.421
##
## P value adjustment method: fdr
pairwise.wilcox.test(KEGG_carbohydrate3$`PTS system, galactosamine-specific IID component`, KEGG_carbohydrate3$Zones,p.adjust.method = "fdr")##
## Pairwise comparisons using Wilcoxon rank sum test
##
## data: KEGG_carbohydrate3$`PTS system, galactosamine-specific IID component` and KEGG_carbohydrate3$Zones
##
## Lam T10
## T10 0.024 -
## T4 0.024 0.690
##
## P value adjustment method: fdr
pairwise.wilcox.test(KEGG_carbohydrate3$`PTS system, lactose-specific IIA component [EC:2.7.1.69]`, KEGG_carbohydrate3$Zones,p.adjust.method = "fdr")##
## Pairwise comparisons using Wilcoxon rank sum test
##
## data: KEGG_carbohydrate3$`PTS system, lactose-specific IIA component [EC:2.7.1.69]` and KEGG_carbohydrate3$Zones
##
## Lam T10
## T10 0.024 -
## T4 0.024 0.548
##
## P value adjustment method: fdr
pairwise.wilcox.test(KEGG_carbohydrate3$`PTS system, N-acetylgalactosamine-specific IIA component`, KEGG_carbohydrate3$Zones,p.adjust.method = "fdr")##
## Pairwise comparisons using Wilcoxon rank sum test
##
## data: KEGG_carbohydrate3$`PTS system, N-acetylgalactosamine-specific IIA component` and KEGG_carbohydrate3$Zones
##
## Lam T10
## T10 0.024 -
## T4 0.024 0.548
##
## P value adjustment method: fdr
pairwise.wilcox.test(KEGG_carbohydrate3$`glucose-1-phosphatase [EC:3.1.3.10]`, KEGG_carbohydrate3$Zones,p.adjust.method = "fdr")##
## Pairwise comparisons using Wilcoxon rank sum test
##
## data: KEGG_carbohydrate3$`glucose-1-phosphatase [EC:3.1.3.10]` and KEGG_carbohydrate3$Zones
##
## Lam T10
## T10 0.024 -
## T4 0.024 0.690
##
## P value adjustment method: fdr
pairwise.wilcox.test(KEGG_carbohydrate3$`phosphoglycerate mutase [EC:5.4.2.1]`, KEGG_carbohydrate3$Zones,p.adjust.method = "fdr")##
## Pairwise comparisons using Wilcoxon rank sum test
##
## data: KEGG_carbohydrate3$`phosphoglycerate mutase [EC:5.4.2.1]` and KEGG_carbohydrate3$Zones
##
## Lam T10
## T10 0.024 -
## T4 0.024 0.421
##
## P value adjustment method: fdr
pairwise.wilcox.test(KEGG_carbohydrate3$`pyruvate dehydrogenase E2 component (dihydrolipoamide`, KEGG_carbohydrate3$Zones,p.adjust.method = "fdr")##
## Pairwise comparisons using Wilcoxon rank sum test
##
## data: KEGG_carbohydrate3$`pyruvate dehydrogenase E2 component (dihydrolipoamide` and KEGG_carbohydrate3$Zones
##
## Lam T10
## T10 0.690 -
## T4 0.024 0.083
##
## P value adjustment method: fdr
pairwise.wilcox.test(KEGG_carbohydrate3$`phosphoenolpyruvate carboxykinase (GTP) [EC:4.1.1.32]`, KEGG_carbohydrate3$Zones,p.adjust.method = "fdr")##
## Pairwise comparisons using Wilcoxon rank sum test
##
## data: KEGG_carbohydrate3$`phosphoenolpyruvate carboxykinase (GTP) [EC:4.1.1.32]` and KEGG_carbohydrate3$Zones
##
## Lam T10
## T10 0.333 -
## T4 0.095 0.841
##
## P value adjustment method: fdr
pairwise.wilcox.test(KEGG_carbohydrate3$`phosphoglucomutase [EC:5.4.2.2]`, KEGG_carbohydrate3$Zones,p.adjust.method = "fdr")##
## Pairwise comparisons using Wilcoxon rank sum test
##
## data: KEGG_carbohydrate3$`phosphoglucomutase [EC:5.4.2.2]` and KEGG_carbohydrate3$Zones
##
## Lam T10
## T10 0.690 -
## T4 0.024 0.690
##
## P value adjustment method: fdr
pairwise.wilcox.test(KEGG_carbohydrate3$`acetaldehyde dehydrogenase / alcohol dehydrogenase [EC:1.2.1.10`, KEGG_carbohydrate3$Zones,p.adjust.method = "fdr")##
## Pairwise comparisons using Wilcoxon rank sum test
##
## data: KEGG_carbohydrate3$`acetaldehyde dehydrogenase / alcohol dehydrogenase [EC:1.2.1.10` and KEGG_carbohydrate3$Zones
##
## Lam T10
## T10 0.333 -
## T4 0.548 0.024
##
## P value adjustment method: fdr
pairwise.wilcox.test(KEGG_carbohydrate3$`S-(hydroxymethyl)glutathione dehydrogenase / alcohol dehydrogenase`, KEGG_carbohydrate3$Zones,p.adjust.method = "fdr")##
## Pairwise comparisons using Wilcoxon rank sum test
##
## data: KEGG_carbohydrate3$`S-(hydroxymethyl)glutathione dehydrogenase / alcohol dehydrogenase` and KEGG_carbohydrate3$Zones
##
## Lam T10
## T10 0.032 -
## T4 0.012 0.012
##
## P value adjustment method: fdr
pairwise.wilcox.test(KEGG_carbohydrate3$`alcohol dehydrogenase [EC:1.1.1.1]`, KEGG_carbohydrate3$Zones,p.adjust.method = "fdr")##
## Pairwise comparisons using Wilcoxon rank sum test
##
## data: KEGG_carbohydrate3$`alcohol dehydrogenase [EC:1.1.1.1]` and KEGG_carbohydrate3$Zones
##
## Lam T10
## T10 0.083 -
## T4 0.083 0.421
##
## P value adjustment method: fdr
pairwise.wilcox.test(KEGG_carbohydrate3$`phosphoglycolate phosphatase [EC:3.1.3.18]`, KEGG_carbohydrate3$Zones,p.adjust.method = "fdr")##
## Pairwise comparisons using Wilcoxon rank sum test
##
## data: KEGG_carbohydrate3$`phosphoglycolate phosphatase [EC:3.1.3.18]` and KEGG_carbohydrate3$Zones
##
## Lam T10
## T10 0.024 -
## T4 0.024 0.548
##
## P value adjustment method: fdr
pairwise.wilcox.test(KEGG_carbohydrate3$`formate dehydrogenase, alpha subunit [EC:1.2.1.2]`, KEGG_carbohydrate3$Zones,p.adjust.method = "fdr")##
## Pairwise comparisons using Wilcoxon rank sum test
##
## data: KEGG_carbohydrate3$`formate dehydrogenase, alpha subunit [EC:1.2.1.2]` and KEGG_carbohydrate3$Zones
##
## Lam T10
## T10 0.690 -
## T4 0.095 0.333
##
## P value adjustment method: fdr
pairwise.wilcox.test(KEGG_carbohydrate3$`formyltetrahydrofolate deformylase [EC:3.5.1.10]`, KEGG_carbohydrate3$Zones,p.adjust.method = "fdr")##
## Pairwise comparisons using Wilcoxon rank sum test
##
## data: KEGG_carbohydrate3$`formyltetrahydrofolate deformylase [EC:3.5.1.10]` and KEGG_carbohydrate3$Zones
##
## Lam T10
## T10 0.048 -
## T4 0.048 0.548
##
## P value adjustment method: fdr
pairwise.wilcox.test(KEGG_carbohydrate3$`myo-inositol-1(or 4)-monophosphatase [EC:3.1.3.25]`, KEGG_carbohydrate3$Zones,p.adjust.method = "fdr")##
## Pairwise comparisons using Wilcoxon rank sum test
##
## data: KEGG_carbohydrate3$`myo-inositol-1(or 4)-monophosphatase [EC:3.1.3.25]` and KEGG_carbohydrate3$Zones
##
## Lam T10
## T10 0.143 -
## T4 0.048 0.841
##
## P value adjustment method: fdr
pairwise.wilcox.test(KEGG_carbohydrate3$`rhamnulose-1-phosphate aldolase [EC:4.1.2.19]`, KEGG_carbohydrate3$Zones,p.adjust.method = "fdr")##
## Pairwise comparisons using Wilcoxon rank sum test
##
## data: KEGG_carbohydrate3$`rhamnulose-1-phosphate aldolase [EC:4.1.2.19]` and KEGG_carbohydrate3$Zones
##
## Lam T10
## T10 0.083 -
## T4 0.083 0.841
##
## P value adjustment method: fdr
pairwise.wilcox.test(KEGG_carbohydrate3$`3-hexulose-6-phosphate synthase [EC:4.1.2.43]`, KEGG_carbohydrate3$Zones,p.adjust.method = "fdr")##
## Pairwise comparisons using Wilcoxon rank sum test
##
## data: KEGG_carbohydrate3$`3-hexulose-6-phosphate synthase [EC:4.1.2.43]` and KEGG_carbohydrate3$Zones
##
## Lam T10
## T10 0.690 -
## T4 0.024 0.048
##
## P value adjustment method: fdr
pairwise.wilcox.test(KEGG_carbohydrate3$`glucose 1-dehydrogenase [EC:1.1.1.47]`, KEGG_carbohydrate3$Zones,p.adjust.method = "fdr")##
## Pairwise comparisons using Wilcoxon rank sum test
##
## data: KEGG_carbohydrate3$`glucose 1-dehydrogenase [EC:1.1.1.47]` and KEGG_carbohydrate3$Zones
##
## Lam T10
## T10 0.548 -
## T4 0.024 0.048
##
## P value adjustment method: fdr
pairwise.wilcox.test(KEGG_carbohydrate3$`ribokinase [EC:2.7.1.15]`, KEGG_carbohydrate3$Zones,p.adjust.method = "fdr")##
## Pairwise comparisons using Wilcoxon rank sum test
##
## data: KEGG_carbohydrate3$`ribokinase [EC:2.7.1.15]` and KEGG_carbohydrate3$Zones
##
## Lam T10
## T10 0.024 -
## T4 0.024 0.841
##
## P value adjustment method: fdr
pairwise.wilcox.test(KEGG_carbohydrate3$`phosphopentomutase [EC:5.4.2.7]`, KEGG_carbohydrate3$Zones,p.adjust.method = "fdr")##
## Pairwise comparisons using Wilcoxon rank sum test
##
## data: KEGG_carbohydrate3$`phosphopentomutase [EC:5.4.2.7]` and KEGG_carbohydrate3$Zones
##
## Lam T10
## T10 0.421 -
## T4 0.421 0.095
##
## P value adjustment method: fdr
pairwise.wilcox.test(KEGG_carbohydrate3$`methylglyoxal synthase [EC:4.2.3.3]`, KEGG_carbohydrate3$Zones,p.adjust.method = "fdr")##
## Pairwise comparisons using Wilcoxon rank sum test
##
## data: KEGG_carbohydrate3$`methylglyoxal synthase [EC:4.2.3.3]` and KEGG_carbohydrate3$Zones
##
## Lam T10
## T10 0.690 -
## T4 0.095 0.631
##
## P value adjustment method: fdr
pairwise.wilcox.test(KEGG_carbohydrate3$`lactaldehyde reductase [EC:1.1.1.77]`, KEGG_carbohydrate3$Zones,p.adjust.method = "fdr")##
## Pairwise comparisons using Wilcoxon rank sum test
##
## data: KEGG_carbohydrate3$`lactaldehyde reductase [EC:1.1.1.77]` and KEGG_carbohydrate3$Zones
##
## Lam T10
## T10 0.226 -
## T4 0.095 1.000
##
## P value adjustment method: fdr
pairwise.wilcox.test(KEGG_carbohydrate3$`acetyl-CoA carboxylase biotin carboxyl carrier protein`, KEGG_carbohydrate3$Zones,p.adjust.method = "fdr")##
## Pairwise comparisons using Wilcoxon rank sum test
##
## data: KEGG_carbohydrate3$`acetyl-CoA carboxylase biotin carboxyl carrier protein` and KEGG_carbohydrate3$Zones
##
## Lam T10
## T10 0.222 -
## T4 0.083 0.024
##
## P value adjustment method: fdr
pairwise.wilcox.test(KEGG_carbohydrate3$`propionate CoA-transferase [EC:2.8.3.1]`, KEGG_carbohydrate3$Zones,p.adjust.method = "fdr")##
## Pairwise comparisons using Wilcoxon rank sum test
##
## data: KEGG_carbohydrate3$`propionate CoA-transferase [EC:2.8.3.1]` and KEGG_carbohydrate3$Zones
##
## Lam T10
## T10 0.548 -
## T4 0.095 0.226
##
## P value adjustment method: fdr
pairwise.wilcox.test(KEGG_carbohydrate3$`phosphoenolpyruvate carboxylase [EC:4.1.1.31]`, KEGG_carbohydrate3$Zones,p.adjust.method = "fdr")##
## Pairwise comparisons using Wilcoxon rank sum test
##
## data: KEGG_carbohydrate3$`phosphoenolpyruvate carboxylase [EC:4.1.1.31]` and KEGG_carbohydrate3$Zones
##
## Lam T10
## T10 0.024 -
## T4 0.024 0.421
##
## P value adjustment method: fdr
pairwise.wilcox.test(KEGG_carbohydrate3$`cellobiose phosphorylase [EC:2.4.1.20]`, KEGG_carbohydrate3$Zones,p.adjust.method = "fdr")##
## Pairwise comparisons using Wilcoxon rank sum test
##
## data: KEGG_carbohydrate3$`cellobiose phosphorylase [EC:2.4.1.20]` and KEGG_carbohydrate3$Zones
##
## Lam T10
## T10 0.048 -
## T4 0.024 0.421
##
## P value adjustment method: fdr
pairwise.wilcox.test(KEGG_carbohydrate3$`maltose phosphorylase [EC:2.4.1.8]`, KEGG_carbohydrate3$Zones,p.adjust.method = "fdr")##
## Pairwise comparisons using Wilcoxon rank sum test
##
## data: KEGG_carbohydrate3$`maltose phosphorylase [EC:2.4.1.8]` and KEGG_carbohydrate3$Zones
##
## Lam T10
## T10 0.048 -
## T4 0.048 0.095
##
## P value adjustment method: fdr
pairwise.wilcox.test(KEGG_carbohydrate3$`arabinan endo-1,5-alpha-L-arabinosidase [EC:3.2.1.99]`, KEGG_carbohydrate3$Zones,p.adjust.method = "fdr")##
## Pairwise comparisons using Wilcoxon rank sum test
##
## data: KEGG_carbohydrate3$`arabinan endo-1,5-alpha-L-arabinosidase [EC:3.2.1.99]` and KEGG_carbohydrate3$Zones
##
## Lam T10
## T10 0.083 -
## T4 0.048 0.548
##
## P value adjustment method: fdr
pairwise.wilcox.test(KEGG_carbohydrate3$`endo-1,4-beta-xylanase [EC:3.2.1.8]`, KEGG_carbohydrate3$Zones,p.adjust.method = "fdr")##
## Pairwise comparisons using Wilcoxon rank sum test
##
## data: KEGG_carbohydrate3$`endo-1,4-beta-xylanase [EC:3.2.1.8]` and KEGG_carbohydrate3$Zones
##
## Lam T10
## T10 0.143 -
## T4 0.048 0.841
##
## P value adjustment method: fdr
pairwise.wilcox.test(KEGG_carbohydrate3$`myo-inositol catabolism protein IolS [EC:1.1.1.-]`, KEGG_carbohydrate3$Zones,p.adjust.method = "fdr")## Warning in wilcox.test.default(xi, xj, paired = paired, ...): cannot
## compute exact p-value with ties
## Warning in wilcox.test.default(xi, xj, paired = paired, ...): cannot
## compute exact p-value with ties
## Warning in wilcox.test.default(xi, xj, paired = paired, ...): cannot
## compute exact p-value with ties
##
## Pairwise comparisons using Wilcoxon rank sum test
##
## data: KEGG_carbohydrate3$`myo-inositol catabolism protein IolS [EC:1.1.1.-]` and KEGG_carbohydrate3$Zones
##
## Lam T10
## T10 0.424 -
## T4 0.015 0.015
##
## P value adjustment method: fdr
pairwise.wilcox.test(KEGG_carbohydrate3$`putative family 31 glucosidase`, KEGG_carbohydrate3$Zones,p.adjust.method = "fdr")##
## Pairwise comparisons using Wilcoxon rank sum test
##
## data: KEGG_carbohydrate3$`putative family 31 glucosidase` and KEGG_carbohydrate3$Zones
##
## Lam T10
## T10 0.143 -
## T4 0.048 0.690
##
## P value adjustment method: fdr
pairwise.wilcox.test(KEGG_carbohydrate3$`sugar fermentation stimulation protein A`, KEGG_carbohydrate3$Zones,p.adjust.method = "fdr")##
## Pairwise comparisons using Wilcoxon rank sum test
##
## data: KEGG_carbohydrate3$`sugar fermentation stimulation protein A` and KEGG_carbohydrate3$Zones
##
## Lam T10
## T10 0.548 -
## T4 0.095 0.333
##
## P value adjustment method: fdr
KEGG_energy=read.csv("mean_KEGGFUN_wilcox_0.05_energy_metabolism.csv",header=T,sep=",")
rownames(KEGG_energy)=KEGG_energy[,1]
KEGG_energy1=KEGG_energy[,-(1:2)]
KEGG_energy1_norm <- t(apply(KEGG_energy1, 1, cal_z_score))
my_sample_col <- data.frame(sample = rep(c("Lam","T10","T4"), each=5))
row.names(my_sample_col) <- colnames(KEGG_energy1_norm)
pheatmap(KEGG_energy1_norm,annotation_col=my_sample_col,show_colnames = F, show_rownames = T,legend = TRUE, cluster_rows = T,cluster_cols = F,border_color = NA,cellwidth = 12, cellheight = 15, main="Energy Metabolism")#Lam and T4
KEGG_energy_functions1 <- KEGG_energy1[,c(idx_Lam, idx_T4)]
clin.sub_KEGG_energy_functions1 <- data.frame(ID = 1:ncol(KEGG_energy_functions1), group = rep(c("Lam","T4"), each=5))
fx <- function(x) {
res <-wilcox.test(x ~clin.sub_KEGG_energy_functions1$group)
return(res$p.value)
}
KEGG_energy_functions.wilcox1 <-apply(as.matrix(KEGG_energy_functions1), 1, fx)
KEGG_energy_functions.wilcox1 <- data.frame(KEGG_energy_functions.wilcox1 )
colnames(KEGG_energy_functions.wilcox1) <- c("p.value")
write.csv(KEGG_energy_functions.wilcox1, file = "KEGG_energy_functions.wilcox0.05_Lam_T4.csv")
#Lam and T10
KEGG_energy_functions2 <- KEGG_energy1[,c(idx_Lam, idx_T10)]
clin.sub_KEGG_energy_functions2 <- data.frame(ID = 1:ncol(KEGG_energy_functions2), group = rep(c("Lam","T10"), each=5))
fx <- function(x) {
res <-wilcox.test(x ~clin.sub_KEGG_energy_functions2$group)
return(res$p.value)
}
KEGG_energy_functions.wilcox2 <-apply(as.matrix(KEGG_energy_functions2), 1, fx)
KEGG_energy_functions.wilcox2 <- data.frame(KEGG_energy_functions.wilcox2 )
colnames(KEGG_energy_functions.wilcox2) <- c("p.value")
write.csv(KEGG_energy_functions.wilcox2, file = "KEGG_energy_functions.wilcox0.05_Lam_T10.csv")
#T10 and T4
KEGG_energy_functions3 <-KEGG_energy1[,c(idx_T4, idx_T10)]
clin.sub_KEGG_energy_functions3 <- data.frame(ID = 1:ncol(KEGG_energy_functions3), group = rep(c("T4","T10"), each=5))
fx <- function(x) {
res <-wilcox.test(x ~clin.sub_KEGG_energy_functions3$group)
return(res$p.value)
}
KEGG_energy_functions.wilcox3 <-apply(as.matrix(KEGG_energy_functions3), 1, fx)
KEGG_energy_functions.wilcox3 <- data.frame(KEGG_energy_functions.wilcox3 )
colnames(KEGG_energy_functions.wilcox3) <- c("p.value")
write.csv(KEGG_energy_functions.wilcox3, file = "KEGG_energy_functions.wilcox0.05_T4_T10.csv")KEGG_energy2=t(KEGG_energy1)
KEGG_energy3=cbind(KEGG_energy2, Colvec_expedition)
pairwise.wilcox.test(KEGG_energy3$`methylenetetrahydrofolate reductase (NADPH) [EC:1.5.1.20]`, KEGG_energy3$Zones,p.adjust.method = "fdr")##
## Pairwise comparisons using Wilcoxon rank sum test
##
## data: KEGG_energy3$`methylenetetrahydrofolate reductase (NADPH) [EC:1.5.1.20]` and KEGG_energy3$Zones
##
## Lam T10
## T10 0.226 -
## T4 0.095 0.841
##
## P value adjustment method: fdr
pairwise.wilcox.test(KEGG_energy3$`carbonic anhydrase [EC:4.2.1.1]`, KEGG_energy3$Zones,p.adjust.method = "fdr")##
## Pairwise comparisons using Wilcoxon rank sum test
##
## data: KEGG_energy3$`carbonic anhydrase [EC:4.2.1.1]` and KEGG_energy3$Zones
##
## Lam T10
## T10 0.083 -
## T4 0.048 0.690
##
## P value adjustment method: fdr
pairwise.wilcox.test(KEGG_energy3$`Nif-specific regulatory protein`, KEGG_energy3$Zones,p.adjust.method = "fdr")##
## Pairwise comparisons using Wilcoxon rank sum test
##
## data: KEGG_energy3$`Nif-specific regulatory protein` and KEGG_energy3$Zones
##
## Lam T10
## T10 0.024 -
## T4 0.222 0.048
##
## P value adjustment method: fdr
pairwise.wilcox.test(KEGG_energy3$`glutamate synthase (NADPH/NADH) large chain [EC:1.4.1.13 1.4.1.14]`, KEGG_energy3$Zones,p.adjust.method = "fdr")##
## Pairwise comparisons using Wilcoxon rank sum test
##
## data: KEGG_energy3$`glutamate synthase (NADPH/NADH) large chain [EC:1.4.1.13 1.4.1.14]` and KEGG_energy3$Zones
##
## Lam T10
## T10 0.024 -
## T4 0.024 0.690
##
## P value adjustment method: fdr
pairwise.wilcox.test(KEGG_energy3$`aspartate--ammonia ligase [EC:6.3.1.1]`, KEGG_energy3$Zones,p.adjust.method = "fdr")##
## Pairwise comparisons using Wilcoxon rank sum test
##
## data: KEGG_energy3$`aspartate--ammonia ligase [EC:6.3.1.1]` and KEGG_energy3$Zones
##
## Lam T10
## T10 0.048 -
## T4 0.024 0.310
##
## P value adjustment method: fdr
pairwise.wilcox.test(KEGG_energy3$`cystathionine beta-lyase [EC:4.4.1.8]`, KEGG_energy3$Zones,p.adjust.method = "fdr")##
## Pairwise comparisons using Wilcoxon rank sum test
##
## data: KEGG_energy3$`cystathionine beta-lyase [EC:4.4.1.8]` and KEGG_energy3$Zones
##
## Lam T10
## T10 0.310 -
## T4 0.048 0.048
##
## P value adjustment method: fdr
pairwise.wilcox.test(KEGG_energy3$`NADH dehydrogenase I subunit F [EC:1.6.5.3]`, KEGG_energy3$Zones,p.adjust.method = "fdr")##
## Pairwise comparisons using Wilcoxon rank sum test
##
## data: KEGG_energy3$`NADH dehydrogenase I subunit F [EC:1.6.5.3]` and KEGG_energy3$Zones
##
## Lam T10
## T10 0.548 -
## T4 0.143 0.048
##
## P value adjustment method: fdr
pairwise.wilcox.test(KEGG_energy3$`V-type H+-transporting ATPase subunit A [EC:3.6.3.14]`, KEGG_energy3$Zones,p.adjust.method = "fdr")##
## Pairwise comparisons using Wilcoxon rank sum test
##
## data: KEGG_energy3$`V-type H+-transporting ATPase subunit A [EC:3.6.3.14]` and KEGG_energy3$Zones
##
## Lam T10
## T10 0.310 -
## T4 0.048 0.310
##
## P value adjustment method: fdr
pairwise.wilcox.test(KEGG_energy3$`V-type H+-transporting ATPase subunit E [EC:3.6.3.14]`, KEGG_energy3$Zones,p.adjust.method = "fdr")##
## Pairwise comparisons using Wilcoxon rank sum test
##
## data: KEGG_energy3$`V-type H+-transporting ATPase subunit E [EC:3.6.3.14]` and KEGG_energy3$Zones
##
## Lam T10
## T10 0.548 -
## T4 0.095 0.333
##
## P value adjustment method: fdr
pairwise.wilcox.test(KEGG_energy3$`F-type H+-transporting ATPase subunit c [EC:3.6.3.14]`, KEGG_energy3$Zones,p.adjust.method = "fdr")##
## Pairwise comparisons using Wilcoxon rank sum test
##
## data: KEGG_energy3$`F-type H+-transporting ATPase subunit c [EC:3.6.3.14]` and KEGG_energy3$Zones
##
## Lam T10
## T10 0.143 -
## T4 0.095 0.548
##
## P value adjustment method: fdr
pairwise.wilcox.test(KEGG_energy3$`cystathionine gamma-synthase [EC:2.5.1.48]`, KEGG_energy3$Zones,p.adjust.method = "fdr")##
## Pairwise comparisons using Wilcoxon rank sum test
##
## data: KEGG_energy3$`cystathionine gamma-synthase [EC:2.5.1.48]` and KEGG_energy3$Zones
##
## Lam T10
## T10 0.095 -
## T4 0.083 0.048
##
## P value adjustment method: fdr
pairwise.wilcox.test(KEGG_energy3$`adenylylsulfate kinase [EC:2.7.1.25]`, KEGG_energy3$Zones,p.adjust.method = "fdr")##
## Pairwise comparisons using Wilcoxon rank sum test
##
## data: KEGG_energy3$`adenylylsulfate kinase [EC:2.7.1.25]` and KEGG_energy3$Zones
##
## Lam T10
## T10 0.421 -
## T4 0.095 0.421
##
## P value adjustment method: fdr
pairwise.wilcox.test(KEGG_energy3$`dihydroorotate dehydrogenase electron transfer subunit`, KEGG_energy3$Zones,p.adjust.method = "fdr")##
## Pairwise comparisons using Wilcoxon rank sum test
##
## data: KEGG_energy3$`dihydroorotate dehydrogenase electron transfer subunit` and KEGG_energy3$Zones
##
## Lam T10
## T10 0.222 -
## T4 0.024 0.222
##
## P value adjustment method: fdr
pairwise.wilcox.test(KEGG_energy3$`electron transfer flavoprotein alpha subunit`, KEGG_energy3$Zones,p.adjust.method = "fdr")##
## Pairwise comparisons using Wilcoxon rank sum test
##
## data: KEGG_energy3$`electron transfer flavoprotein alpha subunit` and KEGG_energy3$Zones
##
## Lam T10
## T10 0.048 -
## T4 0.024 0.690
##
## P value adjustment method: fdr
pairwise.wilcox.test(KEGG_energy3$`electron transport complex protein RnfG`, KEGG_energy3$Zones,p.adjust.method = "fdr")##
## Pairwise comparisons using Wilcoxon rank sum test
##
## data: KEGG_energy3$`electron transport complex protein RnfG` and KEGG_energy3$Zones
##
## Lam T10
## T10 0.012 -
## T4 0.012 0.690
##
## P value adjustment method: fdr
pairwise.wilcox.test(KEGG_energy3$`flavodoxin I`, KEGG_energy3$Zones,p.adjust.method = "fdr")##
## Pairwise comparisons using Wilcoxon rank sum test
##
## data: KEGG_energy3$`flavodoxin I` and KEGG_energy3$Zones
##
## Lam T10
## T10 0.083 -
## T4 0.048 0.548
##
## P value adjustment method: fdr
pairwise.wilcox.test(KEGG_energy3$`hydrogenase-4 component B [EC:1.-.-.-]`, KEGG_energy3$Zones,p.adjust.method = "fdr")##
## Pairwise comparisons using Wilcoxon rank sum test
##
## data: KEGG_energy3$`hydrogenase-4 component B [EC:1.-.-.-]` and KEGG_energy3$Zones
##
## Lam T10
## T10 1.000 -
## T4 0.095 0.143
##
## P value adjustment method: fdr
pairwise.wilcox.test(KEGG_energy3$`NADPH2:quinone reductase [EC:1.6.5.5]`, KEGG_energy3$Zones,p.adjust.method = "fdr")##
## Pairwise comparisons using Wilcoxon rank sum test
##
## data: KEGG_energy3$`NADPH2:quinone reductase [EC:1.6.5.5]` and KEGG_energy3$Zones
##
## Lam T10
## T10 0.048 -
## T4 0.048 0.151
##
## P value adjustment method: fdr
pairwise.wilcox.test(KEGG_energy3$`phenylacetic acid degradation protein`, KEGG_energy3$Zones,p.adjust.method = "fdr")##
## Pairwise comparisons using Wilcoxon rank sum test
##
## data: KEGG_energy3$`phenylacetic acid degradation protein` and KEGG_energy3$Zones
##
## Lam T10
## T10 0.083 -
## T4 0.083 0.222
##
## P value adjustment method: fdr
pairwise.wilcox.test(KEGG_energy3$`Cu2+-exporting ATPase [EC:3.6.3.4]`, KEGG_energy3$Zones,p.adjust.method = "fdr")##
## Pairwise comparisons using Wilcoxon rank sum test
##
## data: KEGG_energy3$`Cu2+-exporting ATPase [EC:3.6.3.4]` and KEGG_energy3$Zones
##
## Lam T10
## T10 0.226 -
## T4 0.095 1.000
##
## P value adjustment method: fdr
KEGG_amino_acid=read.csv("mean_KEGGFUN_wilcox_0.05_amino.acids.csv",header=T,sep=",")
rownames(KEGG_amino_acid)=KEGG_amino_acid[,1]
KEGG_amino_acid1=KEGG_amino_acid[,-(1:2)]
KEGG_amino_acid1_norm <- t(apply(KEGG_amino_acid1, 1, cal_z_score))
my_sample_col <- data.frame(sample = rep(c("Lam","T10","T4"), each=5))
row.names(my_sample_col) <- colnames(KEGG_amino_acid1_norm)
pheatmap(KEGG_amino_acid1_norm,annotation_col=my_sample_col,show_colnames = F, show_rownames = T,legend = TRUE, cluster_rows = T,cluster_cols = F,border_color = NA,cellwidth = 12, cellheight = 15,main="Amino Acid Metabolism")#Lam and T4
KEGG_amino_acid_functions1 <- KEGG_amino_acid1[,c(idx_Lam, idx_T4)]
clin.sub_KEGG_amino_acid_functions1 <- data.frame(ID = 1:ncol(KEGG_amino_acid_functions1), group = rep(c("Lam","T4"), each=5))
fx <- function(x) {
res <-wilcox.test(x ~clin.sub_KEGG_amino_acid_functions1$group)
return(res$p.value)
}
KEGG_amino_acid_functions.wilcox1 <-apply(as.matrix(KEGG_amino_acid_functions1), 1, fx)
KEGG_amino_acid_functions.wilcox1 <- data.frame(KEGG_amino_acid_functions.wilcox1 )
colnames(KEGG_amino_acid_functions.wilcox1) <- c("p.value")
write.csv(KEGG_amino_acid_functions.wilcox1, file = "KEGG_amino_acid_functions.wilcox0.05_Lam_T4.csv")
#Lam and T10
KEGG_amino_acid_functions2 <- KEGG_amino_acid1[,c(idx_Lam, idx_T10)]
clin.sub_KEGG_amino_acid_functions2 <- data.frame(ID = 1:ncol(KEGG_amino_acid_functions2), group = rep(c("Lam","T10"), each=5))
fx <- function(x) {
res <-wilcox.test(x ~clin.sub_KEGG_amino_acid_functions2$group)
return(res$p.value)
}
KEGG_amino_acid_functions.wilcox2 <-apply(as.matrix(KEGG_amino_acid_functions2), 1, fx)
KEGG_amino_acid_functions.wilcox2 <- data.frame(KEGG_amino_acid_functions.wilcox2 )
colnames(KEGG_amino_acid_functions.wilcox2) <- c("p.value")
write.csv(KEGG_amino_acid_functions.wilcox2, file = "KEGG_amino_acid_functions.wilcox0.05_Lam_T10.csv")
#T10 and T4
KEGG_amino_acid_functions3 <-KEGG_amino_acid1[,c(idx_T4, idx_T10)]
clin.sub_KEGG_amino_acid_functions3 <- data.frame(ID = 1:ncol(KEGG_amino_acid_functions3), group = rep(c("T4","T10"), each=5))
fx <- function(x) {
res <-wilcox.test(x ~clin.sub_KEGG_amino_acid_functions3$group)
return(res$p.value)
}
KEGG_amino_acid_functions.wilcox3 <-apply(as.matrix(KEGG_amino_acid_functions3), 1, fx)
KEGG_amino_acid_functions.wilcox3 <- data.frame(KEGG_amino_acid_functions.wilcox3 )
colnames(KEGG_amino_acid_functions.wilcox3) <- c("p.value")
write.csv(KEGG_amino_acid_functions.wilcox3, file = "KEGG_amino_acid_functions.wilcox0.05_T4_T10.csv")#Wilcoxon Rank Sum test with false discovery rates
KEGG_amino_acid2=t(KEGG_amino_acid1)
KEGG_amino_acid3=cbind(KEGG_amino_acid2, Colvec_expedition)
pairwise.wilcox.test(KEGG_amino_acid3$`glutamate 5-kinase [EC:2.7.2.11]`, KEGG_amino_acid3$Zones,p.adjust.method = "fdr")##
## Pairwise comparisons using Wilcoxon rank sum test
##
## data: KEGG_amino_acid3$`glutamate 5-kinase [EC:2.7.2.11]` and KEGG_amino_acid3$Zones
##
## Lam T10
## T10 0.048 -
## T4 0.048 0.690
##
## P value adjustment method: fdr
pairwise.wilcox.test(KEGG_amino_acid3$`glutamate N-acetyltransferase / amino-acid N-acetyltransferase`, KEGG_amino_acid3$Zones,p.adjust.method = "fdr")##
## Pairwise comparisons using Wilcoxon rank sum test
##
## data: KEGG_amino_acid3$`glutamate N-acetyltransferase / amino-acid N-acetyltransferase` and KEGG_amino_acid3$Zones
##
## Lam T10
## T10 0.310 -
## T4 0.095 0.226
##
## P value adjustment method: fdr
pairwise.wilcox.test(KEGG_amino_acid3$`arginine deiminase [EC:3.5.3.6]`, KEGG_amino_acid3$Zones,p.adjust.method = "fdr")##
## Pairwise comparisons using Wilcoxon rank sum test
##
## data: KEGG_amino_acid3$`arginine deiminase [EC:3.5.3.6]` and KEGG_amino_acid3$Zones
##
## Lam T10
## T10 0.012 -
## T4 0.012 0.151
##
## P value adjustment method: fdr
pairwise.wilcox.test(KEGG_amino_acid3$`agmatinase [EC:3.5.3.11]`, KEGG_amino_acid3$Zones,p.adjust.method = "fdr")##
## Pairwise comparisons using Wilcoxon rank sum test
##
## data: KEGG_amino_acid3$`agmatinase [EC:3.5.3.11]` and KEGG_amino_acid3$Zones
##
## Lam T10
## T10 0.841 -
## T4 0.048 0.048
##
## P value adjustment method: fdr
pairwise.wilcox.test(KEGG_amino_acid3$`amidase [EC:3.5.1.4]`, KEGG_amino_acid3$Zones,p.adjust.method = "fdr")##
## Pairwise comparisons using Wilcoxon rank sum test
##
## data: KEGG_amino_acid3$`amidase [EC:3.5.1.4]` and KEGG_amino_acid3$Zones
##
## Lam T10
## T10 0.024 -
## T4 0.024 0.222
##
## P value adjustment method: fdr
pairwise.wilcox.test(KEGG_amino_acid3$`ornithine decarboxylase [EC:4.1.1.17]`, KEGG_amino_acid3$Zones,p.adjust.method = "fdr")##
## Pairwise comparisons using Wilcoxon rank sum test
##
## data: KEGG_amino_acid3$`ornithine decarboxylase [EC:4.1.1.17]` and KEGG_amino_acid3$Zones
##
## Lam T10
## T10 0.024 -
## T4 0.024 1.000
##
## P value adjustment method: fdr
pairwise.wilcox.test(KEGG_amino_acid3$`S-adenosylmethionine decarboxylase [EC:4.1.1.50]`, KEGG_amino_acid3$Zones,p.adjust.method = "fdr")##
## Pairwise comparisons using Wilcoxon rank sum test
##
## data: KEGG_amino_acid3$`S-adenosylmethionine decarboxylase [EC:4.1.1.50]` and KEGG_amino_acid3$Zones
##
## Lam T10
## T10 1.000 -
## T4 0.226 0.095
##
## P value adjustment method: fdr
pairwise.wilcox.test(KEGG_amino_acid3$`D-serine dehydratase [EC:4.3.1.18]`, KEGG_amino_acid3$Zones,p.adjust.method = "fdr")##
## Pairwise comparisons using Wilcoxon rank sum test
##
## data: KEGG_amino_acid3$`D-serine dehydratase [EC:4.3.1.18]` and KEGG_amino_acid3$Zones
##
## Lam T10
## T10 0.143 -
## T4 0.690 0.048
##
## P value adjustment method: fdr
pairwise.wilcox.test(KEGG_amino_acid3$`threonine synthase [EC:4.2.3.1]`, KEGG_amino_acid3$Zones,p.adjust.method = "fdr")##
## Pairwise comparisons using Wilcoxon rank sum test
##
## data: KEGG_amino_acid3$`threonine synthase [EC:4.2.3.1]` and KEGG_amino_acid3$Zones
##
## Lam T10
## T10 0.048 -
## T4 0.048 0.690
##
## P value adjustment method: fdr
pairwise.wilcox.test(KEGG_amino_acid3$`cyclase HisF [EC:4.1.3.-]`, KEGG_amino_acid3$Zones,p.adjust.method = "fdr")##
## Pairwise comparisons using Wilcoxon rank sum test
##
## data: KEGG_amino_acid3$`cyclase HisF [EC:4.1.3.-]` and KEGG_amino_acid3$Zones
##
## Lam T10
## T10 0.083 -
## T4 0.048 1.000
##
## P value adjustment method: fdr
pairwise.wilcox.test(KEGG_amino_acid3$`phosphoribosyl-ATP pyrophosphohydrolase / phosphoribosyl-AMP`, KEGG_amino_acid3$Zones,p.adjust.method = "fdr")##
## Pairwise comparisons using Wilcoxon rank sum test
##
## data: KEGG_amino_acid3$`phosphoribosyl-ATP pyrophosphohydrolase / phosphoribosyl-AMP` and KEGG_amino_acid3$Zones
##
## Lam T10
## T10 0.048 -
## T4 0.048 1.000
##
## P value adjustment method: fdr
pairwise.wilcox.test(KEGG_amino_acid3$`imidazoleglycerol-phosphate dehydratase / histidinol-phosphatase`, KEGG_amino_acid3$Zones,p.adjust.method = "fdr")##
## Pairwise comparisons using Wilcoxon rank sum test
##
## data: KEGG_amino_acid3$`imidazoleglycerol-phosphate dehydratase / histidinol-phosphatase` and KEGG_amino_acid3$Zones
##
## Lam T10
## T10 0.333 -
## T4 0.095 0.548
##
## P value adjustment method: fdr
pairwise.wilcox.test(KEGG_amino_acid3$`phosphoribosyl-ATP pyrophosphohydrolase [EC:3.6.1.31]`, KEGG_amino_acid3$Zones,p.adjust.method = "fdr")##
## Pairwise comparisons using Wilcoxon rank sum test
##
## data: KEGG_amino_acid3$`phosphoribosyl-ATP pyrophosphohydrolase [EC:3.6.1.31]` and KEGG_amino_acid3$Zones
##
## Lam T10
## T10 0.333 -
## T4 0.095 0.690
##
## P value adjustment method: fdr
pairwise.wilcox.test(KEGG_amino_acid3$`phosphoribosyl-AMP cyclohydrolase [EC:3.5.4.19]`, KEGG_amino_acid3$Zones,p.adjust.method = "fdr")##
## Pairwise comparisons using Wilcoxon rank sum test
##
## data: KEGG_amino_acid3$`phosphoribosyl-AMP cyclohydrolase [EC:3.5.4.19]` and KEGG_amino_acid3$Zones
##
## Lam T10
## T10 0.333 -
## T4 0.095 1.000
##
## P value adjustment method: fdr
pairwise.wilcox.test(KEGG_amino_acid3$`2-aminoadipate transaminase [EC:2.6.1.-]`, KEGG_amino_acid3$Zones,p.adjust.method = "fdr")##
## Pairwise comparisons using Wilcoxon rank sum test
##
## data: KEGG_amino_acid3$`2-aminoadipate transaminase [EC:2.6.1.-]` and KEGG_amino_acid3$Zones
##
## Lam T10
## T10 0.226 -
## T4 0.048 0.690
##
## P value adjustment method: fdr
pairwise.wilcox.test(KEGG_amino_acid3$`D-alanine transaminase [EC:2.6.1.21]`, KEGG_amino_acid3$Zones,p.adjust.method = "fdr")##
## Pairwise comparisons using Wilcoxon rank sum test
##
## data: KEGG_amino_acid3$`D-alanine transaminase [EC:2.6.1.21]` and KEGG_amino_acid3$Zones
##
## Lam T10
## T10 0.690 -
## T4 0.083 0.083
##
## P value adjustment method: fdr
pairwise.wilcox.test(KEGG_amino_acid3$`chorismate synthase [EC:4.2.3.5]`, KEGG_amino_acid3$Zones,p.adjust.method = "fdr")##
## Pairwise comparisons using Wilcoxon rank sum test
##
## data: KEGG_amino_acid3$`chorismate synthase [EC:4.2.3.5]` and KEGG_amino_acid3$Zones
##
## Lam T10
## T10 0.048 -
## T4 0.048 0.222
##
## P value adjustment method: fdr
pairwise.wilcox.test(KEGG_amino_acid3$`indole-3-glycerol phosphate synthase [EC:4.1.1.48]`, KEGG_amino_acid3$Zones,p.adjust.method = "fdr")##
## Pairwise comparisons using Wilcoxon rank sum test
##
## data: KEGG_amino_acid3$`indole-3-glycerol phosphate synthase [EC:4.1.1.48]` and KEGG_amino_acid3$Zones
##
## Lam T10
## T10 0.012 -
## T4 0.012 1.000
##
## P value adjustment method: fdr
pairwise.wilcox.test(KEGG_amino_acid3$`anthranilate synthase component II [EC:4.1.3.27]`, KEGG_amino_acid3$Zones,p.adjust.method = "fdr")##
## Pairwise comparisons using Wilcoxon rank sum test
##
## data: KEGG_amino_acid3$`anthranilate synthase component II [EC:4.1.3.27]` and KEGG_amino_acid3$Zones
##
## Lam T10
## T10 0.690 -
## T4 0.024 0.631
##
## P value adjustment method: fdr
pairwise.wilcox.test(KEGG_amino_acid3$`prephenate dehydratase [EC:4.2.1.51]`, KEGG_amino_acid3$Zones,p.adjust.method = "fdr")##
## Pairwise comparisons using Wilcoxon rank sum test
##
## data: KEGG_amino_acid3$`prephenate dehydratase [EC:4.2.1.51]` and KEGG_amino_acid3$Zones
##
## Lam T10
## T10 0.024 -
## T4 0.024 1.000
##
## P value adjustment method: fdr
pairwise.wilcox.test(KEGG_amino_acid3$`D-citramalate synthase [EC:2.3.1.182]`, KEGG_amino_acid3$Zones,p.adjust.method = "fdr")##
## Pairwise comparisons using Wilcoxon rank sum test
##
## data: KEGG_amino_acid3$`D-citramalate synthase [EC:2.3.1.182]` and KEGG_amino_acid3$Zones
##
## Lam T10
## T10 0.083 -
## T4 0.048 0.548
##
## P value adjustment method: fdr
pairwise.wilcox.test(KEGG_amino_acid3$`3-hydroxyisobutyrate dehydrogenase [EC:1.1.1.31]`, KEGG_amino_acid3$Zones,p.adjust.method = "fdr")##
## Pairwise comparisons using Wilcoxon rank sum test
##
## data: KEGG_amino_acid3$`3-hydroxyisobutyrate dehydrogenase [EC:1.1.1.31]` and KEGG_amino_acid3$Zones
##
## Lam T10
## T10 0.548 -
## T4 0.024 0.143
##
## P value adjustment method: fdr
pairwise.wilcox.test(KEGG_amino_acid3$`leucyl aminopeptidase [EC:3.4.11.1]`, KEGG_amino_acid3$Zones,p.adjust.method = "fdr")##
## Pairwise comparisons using Wilcoxon rank sum test
##
## data: KEGG_amino_acid3$`leucyl aminopeptidase [EC:3.4.11.1]` and KEGG_amino_acid3$Zones
##
## Lam T10
## T10 0.421 -
## T4 0.421 0.048
##
## P value adjustment method: fdr
pairwise.wilcox.test(KEGG_amino_acid3$`phosphinothricin acetyltransferase [EC:2.3.1.183]`, KEGG_amino_acid3$Zones,p.adjust.method = "fdr")##
## Pairwise comparisons using Wilcoxon rank sum test
##
## data: KEGG_amino_acid3$`phosphinothricin acetyltransferase [EC:2.3.1.183]` and KEGG_amino_acid3$Zones
##
## Lam T10
## T10 0.024 -
## T4 0.048 0.690
##
## P value adjustment method: fdr
pairwise.wilcox.test(KEGG_amino_acid3$`glyoxylase I family protein`, KEGG_amino_acid3$Zones,p.adjust.method = "fdr")##
## Pairwise comparisons using Wilcoxon rank sum test
##
## data: KEGG_amino_acid3$`glyoxylase I family protein` and KEGG_amino_acid3$Zones
##
## Lam T10
## T10 0.841 -
## T4 0.143 0.095
##
## P value adjustment method: fdr
##lipid metabolism
KEGG_lipid=read.csv("mean_KEGGFUN_wilcox_0.05_Lipid.Metabolism.csv",header=T,sep=",")
rownames(KEGG_lipid)=KEGG_lipid[,1]
KEGG_lipid1=KEGG_lipid[,-(1:2)]
KEGG_lipid1_norm <- t(apply(KEGG_lipid1, 1, cal_z_score))
my_sample_col <- data.frame(sample = rep(c("Lam","T10","T4"), each=5))
row.names(my_sample_col) <- colnames(KEGG_lipid1_norm)
pheatmap(KEGG_lipid1_norm,annotation_col=my_sample_col,show_colnames = F, show_rownames = T,legend = TRUE,border_color = NA,cluster_cols = F,cellwidth = 12, cellheight = 15,main="Lipid Metabolism")#Lam and T4
KEGG_lipid_functions1 <- KEGG_lipid1[,c(idx_Lam, idx_T4)]
clin.sub_KEGG_lipid_functions1 <- data.frame(ID = 1:ncol(KEGG_lipid_functions1), group = rep(c("Lam","T4"), each=5))
fx <- function(x) {
res <-wilcox.test(x ~clin.sub_KEGG_lipid_functions1$group)
return(res$p.value)
}
KEGG_lipid_functions.wilcox1 <-apply(as.matrix(KEGG_lipid_functions1), 1, fx)
KEGG_lipid_functions.wilcox1 <- data.frame(KEGG_lipid_functions.wilcox1 )
colnames(KEGG_lipid_functions.wilcox1) <- c("p.value")
write.csv(KEGG_lipid_functions.wilcox1, file = "KEGG_lipid_functions.wilcox0.05_Lam_T4.csv")
#Lam and T10
KEGG_lipid_functions2 <- KEGG_lipid1[,c(idx_Lam, idx_T10)]
clin.sub_KEGG_lipid_functions2 <- data.frame(ID = 1:ncol(KEGG_lipid_functions2), group = rep(c("Lam","T10"), each=5))
fx <- function(x) {
res <-wilcox.test(x ~clin.sub_KEGG_lipid_functions2$group)
return(res$p.value)
}
KEGG_lipid_functions.wilcox2 <-apply(as.matrix(KEGG_lipid_functions2), 1, fx)
KEGG_lipid_functions.wilcox2 <- data.frame(KEGG_lipid_functions.wilcox2 )
colnames(KEGG_lipid_functions.wilcox2) <- c("p.value")
write.csv(KEGG_lipid_functions.wilcox2, file = "KEGG_lipid_functions.wilcox0.05_Lam_T10.csv")
#T10 and T4
KEGG_lipid_functions3 <-KEGG_lipid1[,c(idx_T4, idx_T10)]
clin.sub_KEGG_lipid_functions3 <- data.frame(ID = 1:ncol(KEGG_lipid_functions3), group = rep(c("T4","T10"), each=5))
fx <- function(x) {
res <-wilcox.test(x ~clin.sub_KEGG_lipid_functions3$group)
return(res$p.value)
}
KEGG_lipid_functions.wilcox3 <-apply(as.matrix(KEGG_lipid_functions3), 1, fx)
KEGG_lipid_functions.wilcox3 <- data.frame(KEGG_lipid_functions.wilcox3 )
colnames(KEGG_lipid_functions.wilcox3) <- c("p.value")
write.csv(KEGG_lipid_functions.wilcox3, file = "KEGG_lipid_functions.wilcox0.05_T4_T10.csv")#Wilcoxon Rank Sum test with "BH" (also known as "fdr")
KEGG_lipid2=t(KEGG_lipid1)
KEGG_lipid3=cbind(KEGG_lipid2, Colvec_expedition)
pairwise.wilcox.test(KEGG_lipid3$`gamma-glutamyltranspeptidase [EC:2.3.2.2]`, KEGG_lipid3$Zones,p.adjust.method = "fdr")##
## Pairwise comparisons using Wilcoxon rank sum test
##
## data: KEGG_lipid3$`gamma-glutamyltranspeptidase [EC:2.3.2.2]` and KEGG_lipid3$Zones
##
## Lam T10
## T10 0.841 -
## T4 0.024 0.024
##
## P value adjustment method: fdr
pairwise.wilcox.test(KEGG_lipid3$`acetyl/propionyl carboxylase subunit alpha`, KEGG_lipid3$Zones,p.adjust.method = "fdr")##
## Pairwise comparisons using Wilcoxon rank sum test
##
## data: KEGG_lipid3$`acetyl/propionyl carboxylase subunit alpha` and KEGG_lipid3$Zones
##
## Lam T10
## T10 0.421 -
## T4 0.024 0.048
##
## P value adjustment method: fdr
pairwise.wilcox.test(KEGG_lipid3$`3R-hydroxymyristoyl ACP dehydrase [EC:4.2.1.-]`, KEGG_lipid3$Zones,p.adjust.method = "fdr")##
## Pairwise comparisons using Wilcoxon rank sum test
##
## data: KEGG_lipid3$`3R-hydroxymyristoyl ACP dehydrase [EC:4.2.1.-]` and KEGG_lipid3$Zones
##
## Lam T10
## T10 0.548 -
## T4 0.048 0.333
##
## P value adjustment method: fdr
pairwise.wilcox.test(KEGG_lipid3$`enoyl-[acyl-carrier protein] reductase I [EC:1.3.1.9]`, KEGG_lipid3$Zones,p.adjust.method = "fdr")##
## Pairwise comparisons using Wilcoxon rank sum test
##
## data: KEGG_lipid3$`enoyl-[acyl-carrier protein] reductase I [EC:1.3.1.9]` and KEGG_lipid3$Zones
##
## Lam T10
## T10 0.048 -
## T4 0.048 1.000
##
## P value adjustment method: fdr
pairwise.wilcox.test(KEGG_lipid3$`rubredoxin-NAD+ reductase [EC:1.18.1.1]`, KEGG_lipid3$Zones,p.adjust.method = "fdr")##
## Pairwise comparisons using Wilcoxon rank sum test
##
## data: KEGG_lipid3$`rubredoxin-NAD+ reductase [EC:1.18.1.1]` and KEGG_lipid3$Zones
##
## Lam T10
## T10 0.048 -
## T4 0.048 1.000
##
## P value adjustment method: fdr
pairwise.wilcox.test(KEGG_lipid3$`1,3-propanediol dehydrogenase [EC:1.1.1.202]`, KEGG_lipid3$Zones,p.adjust.method = "fdr")##
## Pairwise comparisons using Wilcoxon rank sum test
##
## data: KEGG_lipid3$`1,3-propanediol dehydrogenase [EC:1.1.1.202]` and KEGG_lipid3$Zones
##
## Lam T10
## T10 0.048 -
## T4 0.048 0.548
##
## P value adjustment method: fdr
pairwise.wilcox.test(KEGG_lipid3$`dihydroxyacetone kinase, N-terminal domain [EC:2.7.1.-]`, KEGG_lipid3$Zones,p.adjust.method = "fdr")##
## Pairwise comparisons using Wilcoxon rank sum test
##
## data: KEGG_lipid3$`dihydroxyacetone kinase, N-terminal domain [EC:2.7.1.-]` and KEGG_lipid3$Zones
##
## Lam T10
## T10 0.841 -
## T4 0.083 0.083
##
## P value adjustment method: fdr
pairwise.wilcox.test(KEGG_lipid3$`fatty-acyl-CoA synthase [EC:6.2.1.-]`, KEGG_lipid3$Zones,p.adjust.method = "fdr")##
## Pairwise comparisons using Wilcoxon rank sum test
##
## data: KEGG_lipid3$`fatty-acyl-CoA synthase [EC:6.2.1.-]` and KEGG_lipid3$Zones
##
## Lam T10
## T10 0.024 -
## T4 0.024 0.421
##
## P value adjustment method: fdr
pairwise.wilcox.test(KEGG_lipid3$`choloylglycine hydrolase [EC:3.5.1.24]`, KEGG_lipid3$Zones,p.adjust.method = "fdr")##
## Pairwise comparisons using Wilcoxon rank sum test
##
## data: KEGG_lipid3$`choloylglycine hydrolase [EC:3.5.1.24]` and KEGG_lipid3$Zones
##
## Lam T10
## T10 0.083 -
## T4 0.083 0.690
##
## P value adjustment method: fdr
pairwise.wilcox.test(KEGG_lipid3$`arylsulfatase [EC:3.1.6.1]`, KEGG_lipid3$Zones,p.adjust.method = "fdr")##
## Pairwise comparisons using Wilcoxon rank sum test
##
## data: KEGG_lipid3$`arylsulfatase [EC:3.1.6.1]` and KEGG_lipid3$Zones
##
## Lam T10
## T10 0.024 -
## T4 0.222 0.048
##
## P value adjustment method: fdr
##enzyme metabolism
KEGG_enzyme=read.csv("mean_KEGGFUN_wilcox_0.05_enzyme.families.csv",header=T,sep=",")
rownames(KEGG_enzyme)=KEGG_enzyme[,1]
KEGG_enzyme1=KEGG_enzyme[,-(1:2)]
KEGG_enzyme1_norm <- t(apply(KEGG_enzyme1, 1, cal_z_score))
my_sample_col <- data.frame(sample = rep(c("Lam","T10","T4"), each=5))
row.names(my_sample_col) <- colnames(KEGG_enzyme1_norm)
pheatmap(KEGG_enzyme1_norm,annotation_col=my_sample_col,show_colnames = F, show_rownames = T,legend = TRUE,border_color = NA,cluster_cols = F,cellwidth = 12, cellheight = 15, main="Energy Metabolism")#Lam and T4
KEGG_enzyme_functions1 <- KEGG_enzyme1[,c(idx_Lam, idx_T4)]
clin.sub_KEGG_enzyme_functions1 <- data.frame(ID = 1:ncol(KEGG_enzyme_functions1), group = rep(c("Lam","T4"), each=5))
fx <- function(x) {
res <-wilcox.test(x ~clin.sub_KEGG_enzyme_functions1$group)
return(res$p.value)
}
KEGG_enzyme_functions.wilcox1 <-apply(as.matrix(KEGG_enzyme_functions1), 1, fx)
KEGG_enzyme_functions.wilcox1 <- data.frame(KEGG_enzyme_functions.wilcox1 )
colnames(KEGG_enzyme_functions.wilcox1) <- c("p.value")
write.csv(KEGG_enzyme_functions.wilcox1, file = "KEGG_enzyme_functions.wilcox0.05_Lam_T4.csv")
#Lam and T10
KEGG_enzyme_functions2 <- KEGG_enzyme1[,c(idx_Lam, idx_T10)]
clin.sub_KEGG_enzyme_functions2 <- data.frame(ID = 1:ncol(KEGG_enzyme_functions2), group = rep(c("Lam","T10"), each=5))
fx <- function(x) {
res <-wilcox.test(x ~clin.sub_KEGG_enzyme_functions2$group)
return(res$p.value)
}
KEGG_enzyme_functions.wilcox2 <-apply(as.matrix(KEGG_enzyme_functions2), 1, fx)
KEGG_enzyme_functions.wilcox2 <- data.frame(KEGG_enzyme_functions.wilcox2 )
colnames(KEGG_enzyme_functions.wilcox2) <- c("p.value")
write.csv(KEGG_enzyme_functions.wilcox2, file = "KEGG_enzyme_functions.wilcox0.05_Lam_T10.csv")
#T10 and T4
KEGG_enzyme_functions3 <-KEGG_enzyme1[,c(idx_T4, idx_T10)]
clin.sub_KEGG_enzyme_functions3 <- data.frame(ID = 1:ncol(KEGG_enzyme_functions3), group = rep(c("T4","T10"), each=5))
fx <- function(x) {
res <-wilcox.test(x ~clin.sub_KEGG_enzyme_functions3$group)
return(res$p.value)
}
KEGG_enzyme_functions.wilcox3 <-apply(as.matrix(KEGG_enzyme_functions3), 1, fx)
KEGG_enzyme_functions.wilcox3 <- data.frame(KEGG_enzyme_functions.wilcox3 )
colnames(KEGG_enzyme_functions.wilcox3) <- c("p.value")
write.csv(KEGG_enzyme_functions.wilcox3, file = "KEGG_enzyme_functions.wilcox0.05_T4_T10.csv")KEGG_enzyme2=t(KEGG_enzyme1)
KEGG_enzyme3=cbind(KEGG_enzyme2, Colvec_expedition)
pairwise.wilcox.test(KEGG_enzyme3$`ATP-dependent Lon protease [EC:3.4.21.53]`, KEGG_lipid3$Zones,p.adjust.method = "fdr")##
## Pairwise comparisons using Wilcoxon rank sum test
##
## data: KEGG_enzyme3$`ATP-dependent Lon protease [EC:3.4.21.53]` and KEGG_lipid3$Zones
##
## Lam T10
## T10 0.024 -
## T4 0.048 1.000
##
## P value adjustment method: fdr
pairwise.wilcox.test(KEGG_enzyme3$`dipeptidase A [EC:3.4.-.-]`, KEGG_lipid3$Zones,p.adjust.method = "fdr")##
## Pairwise comparisons using Wilcoxon rank sum test
##
## data: KEGG_enzyme3$`dipeptidase A [EC:3.4.-.-]` and KEGG_lipid3$Zones
##
## Lam T10
## T10 0.048 -
## T4 0.048 0.151
##
## P value adjustment method: fdr
pairwise.wilcox.test(KEGG_enzyme3$`prolyl oligopeptidase [EC:3.4.21.26]`, KEGG_lipid3$Zones,p.adjust.method = "fdr")##
## Pairwise comparisons using Wilcoxon rank sum test
##
## data: KEGG_enzyme3$`prolyl oligopeptidase [EC:3.4.21.26]` and KEGG_lipid3$Zones
##
## Lam T10
## T10 0.048 -
## T4 0.048 0.421
##
## P value adjustment method: fdr
pairwise.wilcox.test(KEGG_enzyme3$`protease I [EC:3.2.-.-]`, KEGG_lipid3$Zones,p.adjust.method = "fdr")##
## Pairwise comparisons using Wilcoxon rank sum test
##
## data: KEGG_enzyme3$`protease I [EC:3.2.-.-]` and KEGG_lipid3$Zones
##
## Lam T10
## T10 0.083 -
## T4 0.048 0.222
##
## P value adjustment method: fdr
pairwise.wilcox.test(KEGG_enzyme3$`pyroglutamyl-peptidase [EC:3.4.19.3]`, KEGG_lipid3$Zones,p.adjust.method = "fdr")##
## Pairwise comparisons using Wilcoxon rank sum test
##
## data: KEGG_enzyme3$`pyroglutamyl-peptidase [EC:3.4.19.3]` and KEGG_lipid3$Zones
##
## Lam T10
## T10 0.083 -
## T4 0.024 0.548
##
## P value adjustment method: fdr
pairwise.wilcox.test(KEGG_enzyme3$`sortase B`, KEGG_lipid3$Zones,p.adjust.method = "fdr")##
## Pairwise comparisons using Wilcoxon rank sum test
##
## data: KEGG_enzyme3$`sortase B` and KEGG_lipid3$Zones
##
## Lam T10
## T10 0.548 -
## T4 0.548 0.095
##
## P value adjustment method: fdr
pairwise.wilcox.test(KEGG_enzyme3$`dipeptidyl-peptidase 4 [EC:3.4.14.5]`, KEGG_lipid3$Zones,p.adjust.method = "fdr")##
## Pairwise comparisons using Wilcoxon rank sum test
##
## data: KEGG_enzyme3$`dipeptidyl-peptidase 4 [EC:3.4.14.5]` and KEGG_lipid3$Zones
##
## Lam T10
## T10 0.083 -
## T4 0.048 0.222
##
## P value adjustment method: fdr
pairwise.wilcox.test(KEGG_enzyme3$`lactocepin [EC:3.4.21.96]`, KEGG_lipid3$Zones,p.adjust.method = "fdr")##
## Pairwise comparisons using Wilcoxon rank sum test
##
## data: KEGG_enzyme3$`lactocepin [EC:3.4.21.96]` and KEGG_lipid3$Zones
##
## Lam T10
## T10 0.024 -
## T4 0.024 0.548
##
## P value adjustment method: fdr
pairwise.wilcox.test(KEGG_enzyme3$`two-component system, OmpR family, sensor kinase [EC:2.7.13.3]`, KEGG_lipid3$Zones,p.adjust.method = "fdr")##
## Pairwise comparisons using Wilcoxon rank sum test
##
## data: KEGG_enzyme3$`two-component system, OmpR family, sensor kinase [EC:2.7.13.3]` and KEGG_lipid3$Zones
##
## Lam T10
## T10 0.024 -
## T4 0.024 0.690
##
## P value adjustment method: fdr
pairwise.wilcox.test(KEGG_enzyme3$`two-component system, AgrA family, sensor histidine kinase AgrC`, KEGG_lipid3$Zones,p.adjust.method = "fdr")##
## Pairwise comparisons using Wilcoxon rank sum test
##
## data: KEGG_enzyme3$`two-component system, AgrA family, sensor histidine kinase AgrC` and KEGG_lipid3$Zones
##
## Lam T10
## T10 0.024 -
## T4 0.024 0.841
##
## P value adjustment method: fdr
pairwise.wilcox.test(KEGG_enzyme3$`two-component system, NtrC family, sensor histidine kinase HydH`, KEGG_lipid3$Zones,p.adjust.method = "fdr")##
## Pairwise comparisons using Wilcoxon rank sum test
##
## data: KEGG_enzyme3$`two-component system, NtrC family, sensor histidine kinase HydH` and KEGG_lipid3$Zones
##
## Lam T10
## T10 0.024 -
## T4 0.024 0.222
##
## P value adjustment method: fdr
pairwise.wilcox.test(KEGG_enzyme3$`serine/threonine-protein kinase WNK1 [EC:2.7.11.1]`, KEGG_lipid3$Zones,p.adjust.method = "fdr")##
## Pairwise comparisons using Wilcoxon rank sum test
##
## data: KEGG_enzyme3$`serine/threonine-protein kinase WNK1 [EC:2.7.11.1]` and KEGG_lipid3$Zones
##
## Lam T10
## T10 0.048 -
## T4 0.048 1.000
##
## P value adjustment method: fdr
##cofactors metabolism
KEGG_cofactors=read.csv("mean_KEGGFUN_wilcox_0.05_Metabolism.of.Cofactors.csv",header=T,sep=",")
rownames(KEGG_cofactors)=KEGG_cofactors[,1]
KEGG_cofactors1=KEGG_cofactors[,-(1:2)]
KEGG_cofactors1_norm <- t(apply(KEGG_cofactors1, 1, cal_z_score))
my_sample_col <- data.frame(sample = rep(c("Lam","T10","T4"), each=5))
row.names(my_sample_col) <- colnames(KEGG_cofactors1_norm)
pheatmap(KEGG_cofactors1_norm,annotation_col=my_sample_col,show_colnames = F, show_rownames = T,legend = TRUE,cluster_cols = F,border_color = NA,cellwidth = 12, cellheight = 15,main="Cofactors Metabolism")#Lam and T4
KEGG_cofactors_functions1 <- KEGG_cofactors1[,c(idx_Lam, idx_T4)]
clin.sub_KEGG_cofactors_functions1 <- data.frame(ID = 1:ncol(KEGG_cofactors_functions1), group = rep(c("Lam","T4"), each=5))
fx <- function(x) {
res <-wilcox.test(x ~clin.sub_KEGG_cofactors_functions1$group)
return(res$p.value)
}
KEGG_cofactors_functions.wilcox1 <-apply(as.matrix(KEGG_cofactors_functions1), 1, fx)## Warning in wilcox.test.default(x = structure(c(0, 0, 0.392, 0, 0), .Names =
## c("Lam_4.x", : cannot compute exact p-value with ties
KEGG_cofactors_functions.wilcox1 <- data.frame(KEGG_cofactors_functions.wilcox1 )
colnames(KEGG_cofactors_functions.wilcox1) <- c("p.value")
write.csv(KEGG_cofactors_functions.wilcox1, file = "KEGG_cofactors_functions.wilcox0.05_Lam_T4.csv")
#Lam and T10
KEGG_cofactors_functions2 <- KEGG_cofactors1[,c(idx_Lam, idx_T10)]
clin.sub_KEGG_cofactors_functions2 <- data.frame(ID = 1:ncol(KEGG_cofactors_functions2), group = rep(c("Lam","T10"), each=5))
fx <- function(x) {
res <-wilcox.test(x ~clin.sub_KEGG_cofactors_functions2$group)
return(res$p.value)
}
KEGG_cofactors_functions.wilcox2 <-apply(as.matrix(KEGG_cofactors_functions2), 1, fx)## Warning in wilcox.test.default(x = structure(c(0, 0, 0.392, 0, 0), .Names =
## c("Lam_4.x", : cannot compute exact p-value with ties
KEGG_cofactors_functions.wilcox2 <- data.frame(KEGG_cofactors_functions.wilcox2 )
colnames(KEGG_cofactors_functions.wilcox2) <- c("p.value")
write.csv(KEGG_cofactors_functions.wilcox2, file = "KEGG_cofactors_functions.wilcox0.05_Lam_T10.csv")
#T10 and T4
KEGG_cofactors_functions3 <-KEGG_cofactors1[,c(idx_T4, idx_T10)]
clin.sub_KEGG_cofactors_functions3 <- data.frame(ID = 1:ncol(KEGG_cofactors_functions3), group = rep(c("T4","T10"), each=5))
fx <- function(x) {
res <-wilcox.test(x ~clin.sub_KEGG_cofactors_functions3$group)
return(res$p.value)
}
KEGG_cofactors_functions.wilcox3 <-apply(as.matrix(KEGG_cofactors_functions3), 1, fx)## Warning in wilcox.test.default(x = structure(c(0.13, 0, 0.224, 0,
## 0), .Names = c("T10_12.x", : cannot compute exact p-value with ties
KEGG_cofactors_functions.wilcox3 <- data.frame(KEGG_cofactors_functions.wilcox3 )
colnames(KEGG_cofactors_functions.wilcox3) <- c("p.value")
write.csv(KEGG_cofactors_functions.wilcox3, file = "KEGG_cofactors_functions.wilcox0.05_T4_T10.csv")KEGG_cofactors2=t(KEGG_cofactors1)
KEGG_cofactors3=cbind(KEGG_cofactors2, Colvec_expedition)
pairwise.wilcox.test(KEGG_cofactors3$`dihydroneopterin aldolase [EC:4.1.2.25]`, KEGG_cofactors3$Zones,p.adjust.method = "fdr")##
## Pairwise comparisons using Wilcoxon rank sum test
##
## data: KEGG_cofactors3$`dihydroneopterin aldolase [EC:4.1.2.25]` and KEGG_cofactors3$Zones
##
## Lam T10
## T10 0.024 -
## T4 0.024 1.000
##
## P value adjustment method: fdr
pairwise.wilcox.test(KEGG_cofactors3$`para-aminobenzoate synthetase component II [EC:2.6.1.85]`, KEGG_cofactors3$Zones,p.adjust.method = "fdr")##
## Pairwise comparisons using Wilcoxon rank sum test
##
## data: KEGG_cofactors3$`para-aminobenzoate synthetase component II [EC:2.6.1.85]` and KEGG_cofactors3$Zones
##
## Lam T10
## T10 0.548 -
## T4 0.048 0.048
##
## P value adjustment method: fdr
pairwise.wilcox.test(KEGG_cofactors3$`dihydrofolate reductase [EC:1.5.1.3]`, KEGG_cofactors3$Zones,p.adjust.method = "fdr")##
## Pairwise comparisons using Wilcoxon rank sum test
##
## data: KEGG_cofactors3$`dihydrofolate reductase [EC:1.5.1.3]` and KEGG_cofactors3$Zones
##
## Lam T10
## T10 0.143 -
## T4 0.024 0.151
##
## P value adjustment method: fdr
pairwise.wilcox.test(KEGG_cofactors3$`lipoate-protein ligase A [EC:2.7.7.63]`, KEGG_cofactors3$Zones,p.adjust.method = "fdr")##
## Pairwise comparisons using Wilcoxon rank sum test
##
## data: KEGG_cofactors3$`lipoate-protein ligase A [EC:2.7.7.63]` and KEGG_cofactors3$Zones
##
## Lam T10
## T10 0.421 -
## T4 0.024 0.143
##
## P value adjustment method: fdr
pairwise.wilcox.test(KEGG_cofactors3$`dephospho-CoA kinase [EC:2.7.1.24]`, KEGG_cofactors3$Zones,p.adjust.method = "fdr")##
## Pairwise comparisons using Wilcoxon rank sum test
##
## data: KEGG_cofactors3$`dephospho-CoA kinase [EC:2.7.1.24]` and KEGG_cofactors3$Zones
##
## Lam T10
## T10 0.012 -
## T4 0.012 0.222
##
## P value adjustment method: fdr
pairwise.wilcox.test(KEGG_cofactors3$`type III pantothenate kinase [EC:2.7.1.33]`, KEGG_cofactors3$Zones,p.adjust.method = "fdr")##
## Pairwise comparisons using Wilcoxon rank sum test
##
## data: KEGG_cofactors3$`type III pantothenate kinase [EC:2.7.1.33]` and KEGG_cofactors3$Zones
##
## Lam T10
## T10 0.024 -
## T4 0.024 0.548
##
## P value adjustment method: fdr
pairwise.wilcox.test(KEGG_cofactors3$`adenosylcobinamide-phosphate synthase CobD [EC:6.3.1.10]`, KEGG_cofactors3$Zones,p.adjust.method = "fdr")##
## Pairwise comparisons using Wilcoxon rank sum test
##
## data: KEGG_cofactors3$`adenosylcobinamide-phosphate synthase CobD [EC:6.3.1.10]` and KEGG_cofactors3$Zones
##
## Lam T10
## T10 0.048 -
## T4 0.024 0.841
##
## P value adjustment method: fdr
pairwise.wilcox.test(KEGG_cofactors3$`anaerobic magnesium-protoporphyrin IX monomethyl ester cyclase`, KEGG_cofactors3$Zones,p.adjust.method = "fdr")##
## Pairwise comparisons using Wilcoxon rank sum test
##
## data: KEGG_cofactors3$`anaerobic magnesium-protoporphyrin IX monomethyl ester cyclase` and KEGG_cofactors3$Zones
##
## Lam T10
## T10 0.048 -
## T4 0.024 0.421
##
## P value adjustment method: fdr
pairwise.wilcox.test(KEGG_cofactors3$`cobalamin biosynthesis protein CbiG`, KEGG_cofactors3$Zones,p.adjust.method = "fdr")##
## Pairwise comparisons using Wilcoxon rank sum test
##
## data: KEGG_cofactors3$`cobalamin biosynthesis protein CbiG` and KEGG_cofactors3$Zones
##
## Lam T10
## T10 0.226 -
## T4 0.095 0.690
##
## P value adjustment method: fdr
pairwise.wilcox.test(KEGG_cofactors3$`acid phosphatase (class A) [EC:3.1.3.2]`, KEGG_cofactors3$Zones,p.adjust.method = "fdr")##
## Pairwise comparisons using Wilcoxon rank sum test
##
## data: KEGG_cofactors3$`acid phosphatase (class A) [EC:3.1.3.2]` and KEGG_cofactors3$Zones
##
## Lam T10
## T10 0.012 -
## T4 0.012 0.690
##
## P value adjustment method: fdr
pairwise.wilcox.test(KEGG_cofactors3$`hydroxyethylthiazole kinase [EC:2.7.1.50]`, KEGG_cofactors3$Zones,p.adjust.method = "fdr")##
## Pairwise comparisons using Wilcoxon rank sum test
##
## data: KEGG_cofactors3$`hydroxyethylthiazole kinase [EC:2.7.1.50]` and KEGG_cofactors3$Zones
##
## Lam T10
## T10 0.690 -
## T4 0.012 0.012
##
## P value adjustment method: fdr
pairwise.wilcox.test(KEGG_cofactors3$`2-succinyl-5-enolpyruvyl-6-hydroxy-3-cyclohexene-1-carboxylate`, KEGG_cofactors3$Zones,p.adjust.method = "fdr")##
## Pairwise comparisons using Wilcoxon rank sum test
##
## data: KEGG_cofactors3$`2-succinyl-5-enolpyruvyl-6-hydroxy-3-cyclohexene-1-carboxylate` and KEGG_cofactors3$Zones
##
## Lam T10
## T10 0.143 -
## T4 0.048 0.310
##
## P value adjustment method: fdr
pairwise.wilcox.test(KEGG_cofactors3$`2-succinyl-6-hydroxy-2,4-cyclohexadiene-1-carboxylate synthase`, KEGG_cofactors3$Zones,p.adjust.method = "fdr")## Warning in wilcox.test.default(xi, xj, paired = paired, ...): cannot
## compute exact p-value with ties
## Warning in wilcox.test.default(xi, xj, paired = paired, ...): cannot
## compute exact p-value with ties
## Warning in wilcox.test.default(xi, xj, paired = paired, ...): cannot
## compute exact p-value with ties
##
## Pairwise comparisons using Wilcoxon rank sum test
##
## data: KEGG_cofactors3$`2-succinyl-6-hydroxy-2,4-cyclohexadiene-1-carboxylate synthase` and KEGG_cofactors3$Zones
##
## Lam T10
## T10 0.8 -
## T4 0.1 0.1
##
## P value adjustment method: fdr
pairwise.wilcox.test(KEGG_cofactors3$`O-succinylbenzoate synthase [EC:4.2.1.113]`, KEGG_cofactors3$Zones,p.adjust.method = "fdr")##
## Pairwise comparisons using Wilcoxon rank sum test
##
## data: KEGG_cofactors3$`O-succinylbenzoate synthase [EC:4.2.1.113]` and KEGG_cofactors3$Zones
##
## Lam T10
## T10 0.012 -
## T4 0.012 0.056
##
## P value adjustment method: fdr
pairwise.wilcox.test(KEGG_cofactors3$`O-succinylbenzoic acid--CoA ligase [EC:6.2.1.26]`, KEGG_cofactors3$Zones,p.adjust.method = "fdr")##
## Pairwise comparisons using Wilcoxon rank sum test
##
## data: KEGG_cofactors3$`O-succinylbenzoic acid--CoA ligase [EC:6.2.1.26]` and KEGG_cofactors3$Zones
##
## Lam T10
## T10 0.143 -
## T4 0.048 0.151
##
## P value adjustment method: fdr
pairwise.wilcox.test(KEGG_cofactors3$`pyridoxine 4-dehydrogenase [EC:1.1.1.65]`, KEGG_cofactors3$Zones,p.adjust.method = "fdr")##
## Pairwise comparisons using Wilcoxon rank sum test
##
## data: KEGG_cofactors3$`pyridoxine 4-dehydrogenase [EC:1.1.1.65]` and KEGG_cofactors3$Zones
##
## Lam T10
## T10 0.083 -
## T4 0.024 0.548
##
## P value adjustment method: fdr
pairwise.wilcox.test(KEGG_cofactors3$`starvation sensing protein RspA`, KEGG_cofactors3$Zones,p.adjust.method = "fdr")##
## Pairwise comparisons using Wilcoxon rank sum test
##
## data: KEGG_cofactors3$`starvation sensing protein RspA` and KEGG_cofactors3$Zones
##
## Lam T10
## T10 0.548 -
## T4 0.024 0.024
##
## P value adjustment method: fdr
##glycan biosynthesis and metabolism
KEGG_Glycan=read.csv("mean_KEGGFUN_wilcox_0.05_ Glycan.Biosynthesis.csv",header=T,sep=",")
rownames(KEGG_Glycan)=KEGG_Glycan[,1]
KEGG_Glycan1=KEGG_Glycan[,-(1:2)]
KEGG_Glycan1_norm <- t(apply(KEGG_Glycan1, 1, cal_z_score))
my_sample_col <- data.frame(sample = rep(c("Lam","T10","T4"), each=5))
row.names(my_sample_col) <- colnames(KEGG_Glycan1_norm)
pheatmap(KEGG_Glycan1_norm,annotation_col=my_sample_col,show_colnames = F, show_rownames = T,legend = TRUE,cluster_cols = F,border_color = NA,cellwidth = 12, cellheight = 15,main="Glycan.Biosynthesis and metabolism")#Lam and T4
KEGG_Glycan_functions1 <- KEGG_Glycan1[,c(idx_Lam, idx_T4)]
clin.sub_KEGG_Glycan_functions1 <- data.frame(ID = 1:ncol(KEGG_Glycan_functions1), group = rep(c("Lam","T4"), each=5))
fx <- function(x) {
res <-wilcox.test(x ~clin.sub_KEGG_Glycan_functions1$group)
return(res$p.value)
}
KEGG_Glycan_functions.wilcox1 <-apply(as.matrix(KEGG_Glycan_functions1), 1, fx)
KEGG_Glycan_functions.wilcox1 <- data.frame(KEGG_Glycan_functions.wilcox1 )
colnames(KEGG_Glycan_functions.wilcox1) <- c("p.value")
write.csv(KEGG_Glycan_functions.wilcox1, file = "KEGG_Glycan_functions.wilcox0.05_Lam_T4.csv")
#Lam and T10
KEGG_Glycan_functions2 <- KEGG_Glycan1[,c(idx_Lam, idx_T10)]
clin.sub_KEGG_Glycan_functions2 <- data.frame(ID = 1:ncol(KEGG_Glycan_functions2), group = rep(c("Lam","T10"), each=5))
fx <- function(x) {
res <-wilcox.test(x ~clin.sub_KEGG_Glycan_functions2$group)
return(res$p.value)
}
KEGG_Glycan_functions.wilcox2 <-apply(as.matrix(KEGG_Glycan_functions2), 1, fx)## Warning in wilcox.test.default(x = structure(c(30.129, 46.092, 49.159,
## 50.459, : cannot compute exact p-value with ties
KEGG_Glycan_functions.wilcox2 <- data.frame(KEGG_Glycan_functions.wilcox2 )
colnames(KEGG_Glycan_functions.wilcox2) <- c("p.value")
write.csv(KEGG_Glycan_functions.wilcox2, file = "KEGG_Glycan_functions.wilcox0.05_Lam_T10.csv")
#T10 and T4
KEGG_Glycan_functions3 <-KEGG_Glycan1[,c(idx_T4, idx_T10)]
clin.sub_KEGG_Glycan_functions3 <- data.frame(ID = 1:ncol(KEGG_Glycan_functions3), group = rep(c("T4","T10"), each=5))
fx <- function(x) {
res <-wilcox.test(x ~clin.sub_KEGG_Glycan_functions3$group)
return(res$p.value)
}
KEGG_Glycan_functions.wilcox3 <-apply(as.matrix(KEGG_Glycan_functions3), 1, fx)## Warning in wilcox.test.default(x = structure(c(25.743, 0, 0, 0, 0), .Names
## = c("T10_12.x", : cannot compute exact p-value with ties
KEGG_Glycan_functions.wilcox3 <- data.frame(KEGG_Glycan_functions.wilcox3 )
colnames(KEGG_Glycan_functions.wilcox3) <- c("p.value")
write.csv(KEGG_Glycan_functions.wilcox3, file = "KEGG_Glycan_functions.wilcox0.05_T4_T10.csv")KEGG_Glycan2=t(KEGG_Glycan1)
KEGG_Glycan3=cbind(KEGG_Glycan2, Colvec_expedition)
pairwise.wilcox.test(KEGG_Glycan3$`poly(glycerol-phosphate) alpha-glucosyltransferase [EC:2.4.1.52]`, KEGG_Glycan3$Zones,p.adjust.method = "fdr")##
## Pairwise comparisons using Wilcoxon rank sum test
##
## data: KEGG_Glycan3$`poly(glycerol-phosphate) alpha-glucosyltransferase [EC:2.4.1.52]` and KEGG_Glycan3$Zones
##
## Lam T10
## T10 0.841 -
## T4 0.083 0.083
##
## P value adjustment method: fdr
pairwise.wilcox.test(KEGG_Glycan3$`UDP-N-acetyl-D-mannosaminuronic acid transferase [EC:2.4.1.-]`, KEGG_Glycan3$Zones,p.adjust.method = "fdr")## Warning in wilcox.test.default(xi, xj, paired = paired, ...): cannot
## compute exact p-value with ties
## Warning in wilcox.test.default(xi, xj, paired = paired, ...): cannot
## compute exact p-value with ties
##
## Pairwise comparisons using Wilcoxon rank sum test
##
## data: KEGG_Glycan3$`UDP-N-acetyl-D-mannosaminuronic acid transferase [EC:2.4.1.-]` and KEGG_Glycan3$Zones
##
## Lam T10
## T10 0.015 -
## T4 0.015 0.265
##
## P value adjustment method: fdr
pairwise.wilcox.test(KEGG_Glycan3$`dolichyl-phosphate-mannose-protein mannosyltransferase`, KEGG_Glycan3$Zones,p.adjust.method = "fdr")##
## Pairwise comparisons using Wilcoxon rank sum test
##
## data: KEGG_Glycan3$`dolichyl-phosphate-mannose-protein mannosyltransferase` and KEGG_Glycan3$Zones
##
## Lam T10
## T10 0.151 -
## T4 0.095 0.151
##
## P value adjustment method: fdr
pairwise.wilcox.test(KEGG_Glycan3$`phosphoheptose isomerase [EC:5.-.-.-]`, KEGG_Glycan3$Zones,p.adjust.method = "fdr")##
## Pairwise comparisons using Wilcoxon rank sum test
##
## data: KEGG_Glycan3$`phosphoheptose isomerase [EC:5.-.-.-]` and KEGG_Glycan3$Zones
##
## Lam T10
## T10 0.333 -
## T4 0.095 1.000
##
## P value adjustment method: fdr
pairwise.wilcox.test(KEGG_Glycan3$`galacturonosyltransferase [EC:2.4.1.-]`, KEGG_Glycan3$Zones,p.adjust.method = "fdr")##
## Pairwise comparisons using Wilcoxon rank sum test
##
## data: KEGG_Glycan3$`galacturonosyltransferase [EC:2.4.1.-]` and KEGG_Glycan3$Zones
##
## Lam T10
## T10 1.000 -
## T4 0.048 0.048
##
## P value adjustment method: fdr
pairwise.wilcox.test(KEGG_Glycan3$`undecaprenyl-phosphate galactose phosphotransferase [EC:2.7.8.6]`, KEGG_Glycan3$Zones,p.adjust.method = "fdr")##
## Pairwise comparisons using Wilcoxon rank sum test
##
## data: KEGG_Glycan3$`undecaprenyl-phosphate galactose phosphotransferase [EC:2.7.8.6]` and KEGG_Glycan3$Zones
##
## Lam T10
## T10 0.143 -
## T4 0.095 1.000
##
## P value adjustment method: fdr
pairwise.wilcox.test(KEGG_Glycan3$`alpha-mannosidase [EC:3.2.1.24]`, KEGG_Glycan3$Zones,p.adjust.method = "fdr")##
## Pairwise comparisons using Wilcoxon rank sum test
##
## data: KEGG_Glycan3$`alpha-mannosidase [EC:3.2.1.24]` and KEGG_Glycan3$Zones
##
## Lam T10
## T10 0.226 -
## T4 0.095 0.310
##
## P value adjustment method: fdr
pairwise.wilcox.test(KEGG_Glycan3$`serine/alanine adding enzyme [EC:2.3.2.10]`, KEGG_Glycan3$Zones,p.adjust.method = "fdr")##
## Pairwise comparisons using Wilcoxon rank sum test
##
## data: KEGG_Glycan3$`serine/alanine adding enzyme [EC:2.3.2.10]` and KEGG_Glycan3$Zones
##
## Lam T10
## T10 0.310 -
## T4 0.024 0.024
##
## P value adjustment method: fdr
pairwise.wilcox.test(KEGG_Glycan3$`(heptosyl)LPS beta-1,4-glucosyltransferase [EC:2.4.1.-]`, KEGG_Glycan3$Zones,p.adjust.method = "fdr")##
## Pairwise comparisons using Wilcoxon rank sum test
##
## data: KEGG_Glycan3$`(heptosyl)LPS beta-1,4-glucosyltransferase [EC:2.4.1.-]` and KEGG_Glycan3$Zones
##
## Lam T10
## T10 0.841 -
## T4 0.333 0.095
##
## P value adjustment method: fdr
##other metabolisms
KEGG_others=read.csv("mean_KEGGFUN_wilcox_0.05_Metabolism_others.csv",header=T,sep=",")
rownames(KEGG_others)=KEGG_others[,1]
KEGG_others1=KEGG_others[,-(1:2)]
KEGG_others1_norm <- t(apply(KEGG_others1, 1, cal_z_score))
my_sample_col <- data.frame(sample = rep(c("Lam","T10","T4"), each=5))
row.names(my_sample_col) <- colnames(KEGG_others1_norm)
pheatmap(KEGG_others1_norm,annotation_col=my_sample_col,show_colnames = F, show_rownames = T,legend = TRUE,border_color = NA,cluster_cols = F,cellwidth = 12, cellheight = 15,main="Other Metabolisms")#Lam and T4
KEGG_others_functions1 <- KEGG_others1[,c(idx_Lam, idx_T4)]
clin.sub_KEGG_others_functions1 <- data.frame(ID = 1:ncol(KEGG_others_functions1), group = rep(c("Lam","T4"), each=5))
fx <- function(x) {
res <-wilcox.test(x ~clin.sub_KEGG_others_functions1$group)
return(res$p.value)
}
KEGG_others_functions.wilcox1 <-apply(as.matrix(KEGG_others_functions1), 1, fx)
KEGG_others_functions.wilcox1 <- data.frame(KEGG_others_functions.wilcox1 )
colnames(KEGG_others_functions.wilcox1) <- c("p.value")
write.csv(KEGG_others_functions.wilcox1, file = "KEGG_others_functions.wilcox0.05_Lam_T4.csv")
#Lam and T10
KEGG_others_functions2 <- KEGG_others1[,c(idx_Lam, idx_T10)]
clin.sub_KEGG_others_functions2 <- data.frame(ID = 1:ncol(KEGG_others_functions2), group = rep(c("Lam","T10"), each=5))
fx <- function(x) {
res <-wilcox.test(x ~clin.sub_KEGG_others_functions2$group)
return(res$p.value)
}
KEGG_others_functions.wilcox2 <-apply(as.matrix(KEGG_others_functions2), 1, fx)## Warning in wilcox.test.default(x = structure(c(54.144, 79.048, 84.231,
## 89.995, : cannot compute exact p-value with ties
KEGG_others_functions.wilcox2 <- data.frame(KEGG_others_functions.wilcox2 )
colnames(KEGG_others_functions.wilcox2) <- c("p.value")
write.csv(KEGG_others_functions.wilcox2, file = "KEGG_others_functions.wilcox0.05_Lam_T10.csv")
#T10 and T4
KEGG_others_functions3 <-KEGG_Glycan1[,c(idx_T4, idx_T10)]
clin.sub_KEGG_others_functions3 <- data.frame(ID = 1:ncol(KEGG_others_functions3), group = rep(c("T4","T10"), each=5))
fx <- function(x) {
res <-wilcox.test(x ~clin.sub_KEGG_others_functions3$group)
return(res$p.value)
}
KEGG_others_functions.wilcox3 <-apply(as.matrix(KEGG_others_functions3), 1, fx)## Warning in wilcox.test.default(x = structure(c(25.743, 0, 0, 0, 0), .Names
## = c("T10_12.x", : cannot compute exact p-value with ties
KEGG_others_functions.wilcox3 <- data.frame(KEGG_others_functions.wilcox3 )
colnames(KEGG_others_functions.wilcox3) <- c("p.value")
write.csv(KEGG_others_functions.wilcox3, file = "KEGG_others_functions.wilcox0.05_T4_T10.csv")KEGG_others2=t(KEGG_others1)
KEGG_others3=cbind(KEGG_others2, Colvec_expedition)
pairwise.wilcox.test(KEGG_others3$`2,5-diketo-D-gluconate reductase A [EC:1.1.1.274]`, KEGG_others3$Zones,p.adjust.method = "fdr")##
## Pairwise comparisons using Wilcoxon rank sum test
##
## data: KEGG_others3$`2,5-diketo-D-gluconate reductase A [EC:1.1.1.274]` and KEGG_others3$Zones
##
## Lam T10
## T10 0.024 -
## T4 0.024 1.000
##
## P value adjustment method: fdr
pairwise.wilcox.test(KEGG_others3$`carboxylesterase [EC:3.1.1.1]`, KEGG_others3$Zones,p.adjust.method = "fdr")##
## Pairwise comparisons using Wilcoxon rank sum test
##
## data: KEGG_others3$`carboxylesterase [EC:3.1.1.1]` and KEGG_others3$Zones
##
## Lam T10
## T10 0.048 -
## T4 0.024 0.548
##
## P value adjustment method: fdr
pairwise.wilcox.test(KEGG_others3$`FMN-dependent NADH-azoreductase [EC:1.7.-.-]`, KEGG_others3$Zones,p.adjust.method = "fdr")##
## Pairwise comparisons using Wilcoxon rank sum test
##
## data: KEGG_others3$`FMN-dependent NADH-azoreductase [EC:1.7.-.-]` and KEGG_others3$Zones
##
## Lam T10
## T10 0.095 -
## T4 0.095 0.095
##
## P value adjustment method: fdr
pairwise.wilcox.test(KEGG_others3$`N-acetylglucosaminyldiphosphoundecaprenol [EC:2.4.1.187]`, KEGG_others3$Zones,p.adjust.method = "fdr")##
## Pairwise comparisons using Wilcoxon rank sum test
##
## data: KEGG_others3$`N-acetylglucosaminyldiphosphoundecaprenol [EC:2.4.1.187]` and KEGG_others3$Zones
##
## Lam T10
## T10 0.012 -
## T4 0.012 0.841
##
## P value adjustment method: fdr
pairwise.wilcox.test(KEGG_others3$`putative acetyltransferase [EC:2.3.1.-]`, KEGG_others3$Zones,p.adjust.method = "fdr")##
## Pairwise comparisons using Wilcoxon rank sum test
##
## data: KEGG_others3$`putative acetyltransferase [EC:2.3.1.-]` and KEGG_others3$Zones
##
## Lam T10
## T10 0.310 -
## T4 0.024 0.310
##
## P value adjustment method: fdr
pairwise.wilcox.test(KEGG_others3$`putative colanic acid biosynthesis acetyltransferase WcaF`, KEGG_others3$Zones,p.adjust.method = "fdr")## Warning in wilcox.test.default(xi, xj, paired = paired, ...): cannot
## compute exact p-value with ties
## Warning in wilcox.test.default(xi, xj, paired = paired, ...): cannot
## compute exact p-value with ties
##
## Pairwise comparisons using Wilcoxon rank sum test
##
## data: KEGG_others3$`putative colanic acid biosynthesis acetyltransferase WcaF` and KEGG_others3$Zones
##
## Lam T10
## T10 0.047 -
## T4 0.024 0.265
##
## P value adjustment method: fdr
pairwise.wilcox.test(KEGG_others3$`7-alpha-hydroxysteroid dehydrogenase [EC:1.1.1.159]`, KEGG_others3$Zones,p.adjust.method = "fdr")##
## Pairwise comparisons using Wilcoxon rank sum test
##
## data: KEGG_others3$`7-alpha-hydroxysteroid dehydrogenase [EC:1.1.1.159]` and KEGG_others3$Zones
##
## Lam T10
## T10 0.012 -
## T4 0.012 0.421
##
## P value adjustment method: fdr
pairwise.wilcox.test(KEGG_others3$`phenylacetic acid degradation protein`, KEGG_others3$Zones,p.adjust.method = "fdr")##
## Pairwise comparisons using Wilcoxon rank sum test
##
## data: KEGG_others3$`phenylacetic acid degradation protein` and KEGG_others3$Zones
##
## Lam T10
## T10 0.083 -
## T4 0.083 0.222
##
## P value adjustment method: fdr
pairwise.wilcox.test(KEGG_others3$`4-carboxymuconolactone decarboxylase [EC:4.1.1.44]`, KEGG_others3$Zones,p.adjust.method = "fdr")##
## Pairwise comparisons using Wilcoxon rank sum test
##
## data: KEGG_others3$`4-carboxymuconolactone decarboxylase [EC:4.1.1.44]` and KEGG_others3$Zones
##
## Lam T10
## T10 0.222 -
## T4 0.095 0.222
##
## P value adjustment method: fdr
pairwise.wilcox.test(KEGG_others3$`4-hydroxy-3-methylbut-2-enyl diphosphate reductase [EC:1.17.1.2]`, KEGG_others3$Zones,p.adjust.method = "fdr")##
## Pairwise comparisons using Wilcoxon rank sum test
##
## data: KEGG_others3$`4-hydroxy-3-methylbut-2-enyl diphosphate reductase [EC:1.17.1.2]` and KEGG_others3$Zones
##
## Lam T10
## T10 0.222 -
## T4 0.024 0.222
##
## P value adjustment method: fdr
predicted_pOTUs_genus=read.csv("predicted_hosts_genus_RPKM.csv",header=T,sep=",")
rownames(predicted_pOTUs_genus)=predicted_pOTUs_genus[,1]
predicted_pOTUs_genus1=predicted_pOTUs_genus[,-1]
predicted_pOTUs_genus_norm <- t(apply(predicted_pOTUs_genus1, 1, cal_z_score))
pheatmap(predicted_pOTUs_genus_norm,annotation_col = my_sample_col, show_colnames =T,legend = TRUE,cluster_cols=FALSE)